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Introduction  ond  outline  of  results 


This  report  pursues  two  goals.  These  are 

Cl)  A  comparison  of  currently  used  methods  in  computational  physical 
geodesy.  This  was  the  primary  desire  of  the  contractor. 

(2)  A  feasibility  study  on  the  use  of  the  finite  element  method  for  the 
numerical  solution  of  the  fundamental  problem  of  physical  geodesy.  This 
was  the  primary  desire  of  the  author. 

The  fundamental  problem  of  physical  geodesy  is  the  simultaneous 
determination  of  the  earth's  figure  and  potential  from  geometric  and 
gravimetric  measurements.  The  numerical  solution  requires  a  finite 
parameter i zat i on  of  the  potential  and  -  in  case  of  a  sophisticated 
approach  -  also  of  the  earth's  figure. 

A  comparison  of  various  methods  for  the  detailed  representation  of 
the  earth's  gravity  field  has  recently  been  given  by  Tscherning  (1979). 
His  confirmed  impression  is  that  there  is  a  number  of  competing  methods 
performing  about  equally  well  as  far  as  the  quality  of  results  is 
concerned . 

In  chapter  2  of  the  present  report  various  methods  currently  used 
to  approach  the  problem  of  the  determination  of  the  earths  figure  and 
potential  were  examined  from  the  viewpoint  of  computational  efficiency. 
Methods  like  collocation,  surface  layer,  buried  masspoints,  Bjerhammar's 
method,  lead  to  a  fully  occupied  linear  system  of  equations  to  be 
solved.  The  effort  to  solve  such  a  system  is  proportional  to  N where  N 
is  the  number  of  equations.  Breakdown  due  to  0SU~CPU  times  exceeding  100 
hours  occurs  at  about  N  =  10,000  (this  corresponds  to  a  surface  layer 
solution  with  2° x2°  blocks  near  the  equator) 

If  the  pattern  of  data  and  weights  shows  rotational  symmetry, 
great  savings  in  computation  time  can  be  obtained  by  using  techniques 
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based  on  discrete  Fourier  transform  of  block  circulant  matrices  This 
has  been  shown  by  Colombo  CI98Q).  the  author  feels  that  Colombo's 
approach  is  currently  the  best  one  if  an  essentially  nonr eaundant  set  cf 
surface  data  is  employed.  Such  a  set  is  for  example  given  by  l°xt°  block 
averages  of  gravity  anomalies.  Although  the  quality  of  such  block 
averages  varies  greatly  between  areas,  the  assumption  of  equal  weights 
wil i  not  cause  too  much  harm  to  the  estimated  parameters,  because  there 
is  no  problem  of  adjusting  redundant  data.  The  system  must  take  what  it 
gets  and  hcs  no  choice  to  balance  poor  anomalies  against  better 
observations.  Of  course,  the  accuracy  estimates  obtained  from  such  a 
procedure  are  very  problematic. 

Similar  things  may  be  said  about  GEOFAST  developed  by  TASC  The 
asymptotic  speed  is  even  proportional  to  N  log  N.  The  gain  in  speed  is 
paid  for  by  restricting  applications  to  data  distributed  regularly  on  a 
line  or  within  a  rather  small  plane  rectangle.  Some  possible  trouble 
spots  are  indicated  in  chapter  2.  One  of  them  is  concerned  with 
transporting  a  covariance  from  the  sphere  to  the  plane  Harmonicity  gets 
lost  thereby.  It  would  also  be  interesting  to  have  some  idea  on  the 
propor t i ona I i ty  factor  in  front  of  the  N  log  N  term  estimating  the  CPU 
time. 

Chapters  3  to  6  document  a  feasibility  study  on  the  use  of  the 
finite  element  method  in  physical  geodesy.  This  method  leads  to  a  sparse 
set  of  equations  whose  solution  requires  an  effort  proportional  to  N*, 
where  N  has  the  same  meaning  as  above.  Unfortunately  the  constant  of 
proportionality  is  large.  The  break  even  point  between  the  surface  layer 
and  finite  elements  in  a  global  solution  is  estimated  to  be  around  2°x2° 
blocks.  For  smaller  blocks  finite  elements  are  faster,  for  larger  ones 
the  surface  layer  is  faster.  The  effort  for  a  global  solution  based  on 
l°xl°  gravity  anomaly  data  is  estimated  at  700  OSU  CPU  hours.  A  special 
technique  exploiting  the  remote  zone  effect  could  reduce  this  to  about 
250  hours  (the  surface  layer  method  would  require  15,000  hours).  An 
effort  of  250  OSU  CPU  hours  is  considered  loo  large.  Several  reruns 
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would  be  necessary  before  a  satisfactory  choice  of  weights  is  found 
Although  there  exist  computers,  as  for  example  the  ILLIAC  IV,  on  which 
the  CPU  time  could  be  cut  by  a  factor  of  about  64,  the  problem  appears 
too  large  for  one  individual  researcher  or  a  small  research  group. 

Fortunately,  the  remote  zone  effect  allows  to  compute  local 
solutions.  The  report  will  give  an  estimate  for  the  calculation  of  a 
detailed  potential  in  an  equatorial  strip  C 6 . 5  hours)  and  in  a 
rectangular  area  of  size  32°x64°  (covering  eg.  the  contiguous  US).  In 
the  latter  case  30'x30'  data  were  assumed.  CPU  time  was  estimated  at  20 
OSU  hours.  By  a  sophisticated  use  of  the  remote  zone  effect  this  can 
probably  be  lowered  to  6  hours  This  compares  favorably  with  a  surface 
layer  solution  requiring  about  70  hours. 

The  finite  element  method  does  not  rely  on  any  regular  pattern  of 
observation  and  weights.  (Regularity  could  be  exploited  in  the  same  way 
as  with  the  other  methods.  The  additional  saving  in  CPU  time  would, 
however,  not  be  dramatic.)  In  areas  where  the  field  shows  much  detail, 
smaller  elements  may  be  chosen.  Redundant  data,  as  for  example  gravity 
anomalies  plus  geo i d  heights,  pose  no  problem.  The  method  offers  also 
disadvantages.  Harmonicily  of  the  calculated  field  is  only  approximate. 
The  programming  effort  for  an  efficient  computer  implementation  is 
cons i der ab I e . 

The  finite  element  method  loses  much  of  its  efficiency  if  the  data 
are  not  local .  Local  data  are  composed  of  measurements  taken  in  a  way 
that  one  measurement  involves  only  a  single  point  or  a  small  vicinity  of 
a  point.  A  vicinity  of  a  point  is  considered  small  if  it  contains  only  a 
small  number  of  the  finite  elements  in  its  interior,  A  more  precise 
definition  of  the  locality  of  a  measurement  would  be  that  its 
contribution  to  the  normal  equations  must  not  destroy  the  sparsity 
pattern  resulting  from  a  field  r epr esent at i on  by  means  of  finite 
elements.  Data  obtained  by  integrating  over  an  unknown  orbit  are  not 
local.  Neither  are  misclosures  of  large  inertial  navigation  ioops  There 
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are,  however,  ways  to  incorporate  orbits  at  higher  altitudes  in  an 
ef f i c i ent  way 

The  finite  element  method  lends  itself  to  Helmert  blocking,  or  its 
modern  vcriant,  nested  aissection.  Calculations  for  subregions  (nations, 
continents)  could  be  delegated.  Junction  equations  could  be  combined  at 
a  higher  level,  very  much  in  the  same  way  as  this  is  done  in  continental 
network  adjustment. 

Chapter  6  documents  a  number  of  test  calculations  They  were 
carried  out  with  the  following  goals. 

(1)  To  see  whether  the  Ritz-,  or  the  Trefftz-,  or  the  old  fashioned 
least  squares  principle  should  be  used  (The  latter  is  recommended  for 
the  specific  needs  of  physical  geodesy). 

(2)  To  see  whether  the  use  of  cubic  polynomials  is  sufficient,  or 
whether  quint ic  or  even  higher  degree  polynomials  are  needed  (Cubics 
are  sufficient,  quintics  are  already  hopeless  from  the  CPU  time  point  of 


(3)  To  see  whether  a  certain  type  of  element  partition  making  a  best 
possible  use  of  the  'attenuation  with  altitude  effect"  can  be  employed. 
(The  outcome  was  sat  i  sfactory) 

(4)  To  find  an  appropriate  r epresentat i on  of  the  field  in  the  remote 
outer  space  of  the  earth  (Specially  designed  elements  of  infinite  size 
and  appropriately  chosen  shape  functions  performed  well  in  this 
respect ) . 

(5)  To  get  an  idea  how  the  observational  weights  should  be  balanced 
against  the  weights  odd  I  > ed  to  the  equations  enforcing  the  approximate 
fulfillment  of  Laplace's  equation.  (Reasonable  weights  were  found  by 
experiment.  More  insignt  would  be  desirable) 
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(6)  To  see  whether  a  combination  solution  of  surface  gravity  values  and 
satellite  derived  harmonics  is  possible.  (The  answer  is  yes,  but 
additional  tests  are  necessary  to  identify  procedures  preventing  a 
substantial  increase  in  CPU  time). 

The  experiments  were  carried  out  in  2  dimensions  in  order  to  save  CPU 
time.  Small  scale  3-dimensional  calculations  are  desirable,  but  there 
was  no  time  yet  to  perform  them. 

In  the  authors  opinion  the  finite  element  method  has  a  place  in 
physical  geodesy.  It  is  very  likely  that  a  proposal  by  Junkins  (1979) 
will  be  accepted,  suggesting  to  use  a  finite  element  representat i on  of  a 
completely  known  potential  for  the  purpose  of  rapid  recalculation  in 
real  time  application  and  also  otherwise.  The  authors  feeling  is  that 
finite  elements  are  also  useful  to  parameterise  an  unknown  potential 
during  an  estimation  procedure.  However,  the  method  must  be  cultivated 
somewhat  more  before  a  large  scale  effort  is  attempted.  At  the  end  of 
the  research  period  covered  by  this  report,  the  author  began  to  look 
into  a  hybrid  method  which  combines  finite  elements  with  a  surface  layer 
of  multipoles.  Some  preliminary  statements  on  this  envisioned  method  are 
given  in  chapter  7.  Another  feature  which  makes  finite  elements 
attractive  is  the  possibility  to  attack  in  a  head-on  way  the  free 
boundary  value  problems  of  physical  geodesy.  Some  ideas  how  this  could 
be  accomplished  are  found  in  chapter  5. 
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2.  Review  of  various  methods 

Consider  an  earth-centered  and  earth-fixed  coordinate  system.  Choose  a 
convenient  reference  potential,  e.g  that  one  of  an  equ i potent i a  1 
ellipsoid,  or  one  obtained  from  a  truncated  spherical  harmonics 
expansion.  The  normal  potential  is  assumed  to  completely  absorb  the 
rotational  effect  Hence  the  disturbing  potential  V  is  purely 
gravitational.  In  outer  space  it  satisfies  Laplace's  equation 

AM--  0  (2.1) 

Laplace's  equation  represents  a  local  law.  In  order  to  evaluate  the 
second  order  differential  operator  at  a  certain  location  x,  we  need  only 
information  on  V  in  a  local  neighborhood  of  x. 

The  purpose  of  physical  geodesy  is  the  determination  of  the  earths 
surface  and  potential  from  geodetic  measurements.  There  is  hardly  a  need 
to  point  out  that  the  measurements  are  indirect,  and  that  it  is 
therefore  necessary  to  represent  surface  and  potential  in  terms  of  a 
number  of  unknown  parameters.  Because  computers  perform  only  finitely 
many  operations  in  a  finite  amount  of  time,  the  number  of  parameters 
must  be  finite. 

The  choice  of  an  appropriate  set  of  parameters  is  highly 
nontrivial.  Factors  to  be  taken  into  account  are: 

(#)  Type  of  application,  e.g.,  local  improvements  or  global 
corrections  to  the  field. 

(#)  Type  and  distribution  of  measurements,  in  particular, 
homogeneous  or  heterogeneous  sets  of  data. 

(#)  Mathematical  simplicity  and  elegancy.  Simple  setups  are 
easier  to  program.  Debugging  the  programs  is  less  time 


consum i ng . 


(#)  Computer  time  during  production  runs. 

Let  us  restrict  attention  to  the  potential,  forgetting  temporarily  the 
parameters  describing  the  unknown  reference  surface.  Let  us  review  and 
discuss  representat i ons  of  the  potential  that  have  been  used  by 
geodesists.  Our  main  emphasis  will  be  on  the  last  item  listed  above, 
i.e.  computational  efficiency  during  production  runs. 


2.1.  Spherical  harmonics 

They  satisfy  Laplace's  equation  automatically.  The  trace  functions  with 
respect  to  the  unit  sphere  are  the  surface  spherical  harmonics.  They  are 
the  eigenfunctions  of  any  rotation  invariant  operator  on  the  sphere. 
Confer  Mueller  (1966),  Meissl  (1971a),  Robertson  (1978),  Freeden  (1979) 
for  extensive  discussions.  Spherical  harmonics  provide  great  theoretical 
insight.  They  lead  to  the  concept  of  the  power  spectrum  of  the 
potential.  If  a  problem  can  be  formulated  in  a  way  that  rotation 
symmetry  is  preserved,  then  spherical  harmonics  are  also  of  great 
computational  advantage.  In  this  context  recent  papers  by 
Freeden  (1978),  (1979)  are  pointed  out  where  methods  for  numerical 
integration  of  functions  defined  on  the  sphere  are  specified.  The 
formulas  are  related  to  the  familiar  Gaussian  quadrature  formulas  for 
intervals,  i.e.  they  rely  on  a  weighted  average  of  function  values  at  a 
set  of  discrete  points.  It  is  likely  that  the  formulas  can  be  extended 
to  functions  discretized  in  terms  of  block  averages. 

Spherical  harmonics  are  widely  used  in  satellite  geodesy.  At 
satellite  altitudes  of  1800  km  and  above,  the  potential  is  sufficiently 
attenuated  to  be  properly  represented  by  a  spherical  harmonics  expansion 
of  moderately  large  degree  (N  =  20-30,  or  so). 
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Let  us  now  assume  that  we  are  dealing  with  a  problem  involving  a 
heterogeneous  set  of  data,  suggesting  a  least  squares  setup  by  variation 
of  parameters.  Then  one  oarticular  property  of  spher i ca i  harmonics 
counteracts  computational  efficiency.  This  property  of  spherical 
harmonics  is  that  they  ere  nonzero  almost  everywhere  No  function 
different  from  the  zero  function  can  satisfy  Lap i ace's  equation  and,  at 
the  same  time,  vanish  in  a  part  of  the  domain  having  nonzero  measure 
Thus  spherical  harmonics  fail  to  nave  a  local  support.  If  the  disturbing 
potential  is  represented  in  terms  of  spherical  harmonics  as 
N  *  i 

V(x)»  Z  Z  CijHijCx)  (2.2) 

i :0  j«-l  J  J 

and  if  a  local  measurement  leads  to  a  linear  functional  L(V)  involving 
only  points  in  a  small  vicinity  of  x,  then  we  nevertheless  get  a 
representat i on 

L(V)  *  Z  Z  cu  L (HLj)  (2.3) 

l*o  js-L 

where  most,  if  not  all,  of  the  coefficients  eg  are  nonzero.  Due  to  the 
fa  i  lure  of  the  H ^ j  to  have  local  support,  any  local  measurement  will 
introduce  an  observation  equation  into  the  adjustment  which  has  many 
nonzero  coefficients.  The  normal  equations  will  be  practically  full. 
Solving  a  symmetr i c  fu I  I  system  of  m  equations  (without  complete 
inversion)  requires  about 


(2.4) 


elementary  steps,  one  step  comprising  one  multiplication  followed  by  one 
addition.  Confer  equation  (A. 6)  of  Appendix  A.  Assuming  that  a  computer 
can  perform  about  580,303  such  steps  in  a  second  of  time,  we  arrive  at 
the  following  table,  listing  CPU  times  in  dependence  of  various  choices 
for  m . 


am 


CPU  t i me 


1 00000 


0  3  seconds 
5.6  minutes 
93.0  hours 
10.0  years 


Table  2.1.  CPU  times  for 
solution  of  full  m#m  systems 

m  =  10000  corresponds  to  about  N  =  100  in  the  above  expansion  for  VCx). 
The  actual  time  required  to  solve  very  big  systems  will  be  larger  than 
the  CPU  time  due  to  data  transfer  between  central  and  peripheral  memory. 

It  may  be  a  coincidence  that  present  day  computers  limit  the 
spherical  harmonics  expansion  to  about  N  =  100.  It  also  appears  that  the 
physical  significance  of  sperical  harmonics  coefficients  of  degree 
higher  than  100  is  questionable.  A  local  anomaly  of  the  field  will  have 
a  sperical  harmonics  expansion  with  c^j  tapering  off  to  zero  more 
reluctantly  the  more  pronounced  the  anomaly  is,  i.e.,  the  less  smooth  it 
is.  Hence  a  local  anomaly  is  decomposed  into  components  being  nonlocal. 
This  may  be  desirable  in  such  fields  as  optics  or  acoustics.  In  physical 
geodesy  it  is  undesirable. 


2.2  Surface  layer  representat i on  and  related  methods. 

Under  sufficiently  general  conditions  the  outer  potential  can  be 
represented  by  a  surface  layer.  A  surface  layer  can  be  specified  in 
terms  of  finitely  many  parameters  in  various  ways,  such  as  for  example 
by  constant  values  of  density  in  subregions  composing  the  entire  surface 
of  the  earth.  Confer  Koch  and  Witte  (1971),  Morrison  (1980).  Any  of  the 
surface  elements  generates  a  potential,  i.e.,  a  solution  of  Laplace's 


equation  in  outer  soace .  The  total  potential  is  obtained  by 
superposition.  In  this  respect,  the  presently  discussed  method  does  not 
deviate  from  spherical  harmonics  Also  the  property  of  a  computationally 
rather  unaesirabie  nonlocal  support  carries  over.  Hence  the  above  tcble 
continues  to  give  an  indication  of  the  computational  effort  involved,  if 
m  is  taken  as  the  number  of  surface  elements.  Indeed,  global  solutions 
exceeding  m  =  2003,  which  corresponds  to  about  5  by  5  degree  elements 
near  the  equator,  have  not  been  reported  in  the  literature  (Confer, 
however,  subsection  2.7  below  dealing  with  shortcuts  resulting  from 
symmetrical  configurations). 

The  physical  s i gn i f i cance  of  the  surface  I ayer  is  not  immediate. 
However  they  definitely  offer  the  advantage  of  modelling  local  effects. 
In  areas,  where  the  field  is  very  detailed,  or,  where  a  more  detailed 
know  I edqe  of  the  field  is  desired,  smaller  sized  elements  may  be  chosen. 

Similar  statements  can  be  made  about  buried  masspoints  (confer 
e.g  Needham  (1970)  or  Hardy  and  Goepfert  (1975))  and  about  Bjerhammar's 
method  The  latter  represents  the  potential  by  its  boundary  values  on  a 
sphere  entirely  contained  in  the  interior  of  the  earth.  A  finite 
parameter i zat i on  is  achieved  e.g.  by  par t i t i on i ng  the  surface  of  the 
sphere  into  smail  elements,  and  by  assuming  constant  boundary  values  in 
these  elements.  Confer  Bjerhammar  (1978),  Sjoeberg  (1978). 

The  normal  equation  system  will  be  full,  i.e.  the  above  table 
app . i es .  Local  effects  may  be  conveniently  modelled  by  varying  the 
element  size. 


2.3.  Least  squares  collocation . 
2.3.1.  Krarup's  proposal. 


Although  least  squares  collocation  shares  some  features  with  the  methods 
mentioned  under  2.2,  there  are  some  important  deviations.  Again,  the 


n 


potential  is  represented  as  a  linear  superposition  of  special  solutions 
to  Laplace's  equation: 

n 

v(x)--  H  C;.  L;.K(.><  yL)  C2 .5) 

The  number  of  terms,  however,  now  equals  n,  the  number  of  measurements. 
The  functions 

F-Cv)=  LcK(^yJ  C2.6) 

are  derived  from  a  symmetric  and  positive  definite  kernel  K(x,y)  which 
satisfies  Laplace's  equation  with  respect  to  x.  CDue  to  symmetry,  it 
also  satisfies  Laplace's  equation  with  respect  to  y) .  If  the  i —  t h 
measurement  I i  refers  to  a  functional 

iL  1  LL(V(y))  ;  *  =  *L  (2.7) 

then  (x)  as  given  by  equation  (2.6)  is  taken  as  the  i-th  basis 
funct i on . 

The  function  K(x,y)  is  viewed  as  a  reproducing  kernel.  Thus  it 
defines  a  norm  ||V||.  The  c;.  are  chosen  such  that 

L-t(V(y)W  i  (2.8) 

and  that 

ilV(x)li  *  M in.  (2.9) 

Confer  Krarup  (1969).  The  textbook  by  Moritz  (1980)  may  be  consulted  for 
a  detailed  documentation,  discussion,  presentation  of  extensions,  and 
for  its  bibliography. 

it  is  seen  that  least  squares  col  location  uses  as  many  basis 
functions  as  there  are  observations.  As  compared  to  2.1  and  2.2  this 
number  n  will  frequently  be  larger  than  m,  the  number  of  unknowns  in  the 
other  methods.  Hence  we  can  expect  a  very  good  fit  However,  a  price  has 
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to  be  paid  for  this  The  size  of  the  linear  system  to  be  solved  is  n  by 
n  The  system  is  full  for  the  very  same  reasons  as  given  ear  I  i er :  the 
Fj_(x)'s  satisfy  Laplace's  equation.  Hence  the  computational  effort  is 

nL  (2.10) 

6 

The  mathematical  elegancy  of  ieast  squares  collocation  is  undisputed. 

The  method  eas i ly  takes  any  type  of  heterogeneous  data.  The  choice  of  a 
suitable  covariance  function  is  not  immediate  and  requires  insight. 
Confer  Tschermng  and  Rapp  (1974). 

A  computer  implementation  of  the  collocation  method  is  described 
in  Tschermng  (1978). 

2.3  2.  Least  squares  collocation  using  unknoun  parameters. 

A  somewhat  unsatisfactory  aspect  of  pure  least  squares  collocation  is 
the  following  one.  In  areas  of  insufficient  data  the  predicted  function 
tends  to  approach  the  zero  function.  Since  the  problem  of  physical 
geoaesy  resuits  from  a  linearization  procedure  based  on  the  use  of  a 
reference  surface,  the  consequence  is  that  in  areas  of  insufficient  data 
the  reference  surface  is  predicted.  The  reference  surface  is,  however, 
mostly  chosen  according  to  computational  convenience  rather  than 
according  to  its  approximation  of  the  physical  truth. 

This  unsatisfactory  aspect  can  be  counteracted  by  using  a  setup 
including  unknown  parameters  Confer  Mor i tz  (1980),  section  16.  This 
setup  is  closely  related  to  the  concept  of  generalized  splines  which 
will  be  discussed  farther  below. 
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2.4.  Finite  el ements 

Finite  elements  have  been  very  successful  in  other  disciplines.  There 
exists  an  abundance  of  literature.  As  textbooks  we  mention 
Zienkiewicz  (1971),  Schwarz  (1980),  Ciarlet  (1978),  Strang-Fix  (1973) 
Finite  elements  have  also  been  used  by  some  geodesists.  Cf. 

Szameitat  (1979),  Werner  (1979).  Bosman-Eckhar t-Kub i k  were  early 
geodetic  users  of  finite  element  concepts.  They  applied  piecewise 
polynomials  to  surface  appr ox i mat i on  problems.  The  use  of  finite 
elements  in  Physical  Geodesy  (in  the  narrow  sense)  was  up  to  now 
restricted  to  the  representat i on  of  a  known  potential  for  rapid 
recalculation.  Confer  Junkins  (1977),  (1979),  Engels  (1979).  We  intend 
to  use  the  method  also  during  the  determination  of  the  potential 
together  with  the  earth's  surface.  Our  intended  use  will  be  described  in 
detail  in  the  subsequent  chapters.  In  this  section  we  shall  be  very 
brief. 

The  domain  of  interest  is  subdivided  into  finitely  many 
subregions,  called  elements,  of  preferably  simple  shape.  If  the  region 
is  unbounded,  some  of  the  elements  must  be  of  infinite  size.  We  shall 
mostly  work  with  box-type  elements  partitioning  the  r,cp,  A  parameter 
space  resulting  from  a  choice  of  polar  coordinates.  At  the  boundary  of 
any  element  a  number  of  nodes  is  located.  Any  node  is  shared  by  two  or 
more  elements.  In  our  case  nodes  will  be  mostly  at  the  corners  of  the 
boxes;  but  occasionally  some  are  also  encountered  elsewhere  on  the 
faces . 

To  any  node  a  number  of  parameters  is  associated.  Usually  they 
include  the  value  of  the  potential  there  and  of  some  of  its  derivatives. 
An  interpolation  formula  is  prescribed  which  allows  to  calculate  the 
potential  and  its  derivatives  at  any  point  in  the  interior  or  on  the 
boundary  of  an  element  from  the  parameter  values  of  all  nodes  associated 
with  this  element.  The  interpolating  function  is  analytic  and  of  simple 


shape  m  the  interior  of  the  elements  Across  element  boundaries 
continuity  of  the  function  is  usually  required.  Depending  on  the  type  of 
apolication,  continuity  of  some  derivatives  is  also  needed.  We  propose 
the  use  of  tricubic  polynomials  in  the  elements  such  that  the  resuiting 
potential  is  globally  C*1  continuous  Cit  is  continuous  together  with  its 
first  order  der i vat i ves) . 

If  an  observation  of  tne  potential  refers  to  a  location  within  an 
element,  the  resulting  observation  equation  will  involve  only  parameters 
of  nodes  associated  with  this  element.  As  a  consequence,  the  system  of 
observation  equations  will  be  sparse  and  so  will  be  the  normal  equations 
formed  from  them.  This  is  a  great  computational  advantage.  The  normal 
equations  resulting  from  the  observation  equations  are  not  yet 
suf  f  i  c  lent.  A  potential  represented  by  finite  elements  is  not 
automatically  harmonic.  Harmonicity  must  be  enforced  by  another  set  of 
normal  equations  which  must  be  added  to  the  earlier  ones.  We  shall  call 
this  the  contribution  of  the  field  to  the  normal  equations.  It  is 
obtained  by  minimizing  the  integral  over  the  square  of  the  Laplacean  of 
the  field.  Harmonicity  is  not  fully  ensured  this  way,  but  only 
approximately.  The  sparse  structure  of  the  normal  equations  is  not 
impaired  by  the  field  contribution.  Sparseness  is  the  great  benefit  of 
the  finite  element  method.  If  N  is  the  total  number  of  parameters  in  our 
application  of  the  method,  the  computational  effort  wili  be 

const  N  ^  (2.11) 

As  we  shall  see,  the  total  number  of  parameters  is  appreciably  larger 
than  the  number  of  blocks  in  a  comparable  surface  layer  solution. 

However  it  is  important  to  stress  that  a  synchronized  refinement  of  the 
partition  in  the  two  methods  will  keep  the  ratio  of  the  number  of  blocks 
to  the  number  of  parameters  approximately  constant.  Hence  our  method  is 
constp#NT  as  compared  to  consts *N3  in  the  surface  layer  solution.  On  the 
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other  hand  constF  >  consts  .  As  a  consequence,  for  a  sma I  I  number  of 
blocks  the  surface  layer  solution  will  be  more  economical.  For  large  N, 
the  finite  element  solution  is  better.  The  break  even  point  is  estimated 
to  be  around  2°x2°  blocks. 

Remark •  A  computat i ona I  effort  of  OCN^)  steps  results  if  the  'nested 
dissection  method"  due  to  George  (1973),  (1977)  is  applied  to  a  sparse 
system  of  N  equations  resulting  from  decomposing  a  two  dimensional 
region  of  size  QO/N)xOO/N)  into  N  elements.  The  system  represents  the 
equilibrium  equation  for  a  2-dimensional  elastic  problem  defined  over 
this  region.  Although  our  region  is  3-dimensional,  the  estimate  of  O(N^) 
steps  remains  valid  due  to  the  opportunity  to  use  increasingly  Icrger 
elements  as  the  altitude  increases. 

Remark :  it  must  be  emphasized  that  the  efficiency  of  the  finite  element 
method  relies  heavily  on  the  locality  of  the  measurement  functionals. 

Any  measurement  must  involve  only  one  point  or  a  set  of  points  confined 

to  a  small  region.  Measurements  involving  points  along  an  unknown 
trajectory,  such  as  for  example  misciosures  of  large  inertial  navigation 
loops,  are  excluded.  They  would  destroy  the  sparsity  pattern  and  degrade 
the  asymptotic  computational  efficiency  to  that  of  least  squares 
collocation.  Also  unknown  orbits  of  satellites  pose  difficulties.  If  the 
orbits  are  high  enough,  one  may,  however  compromize  by  fusing  finite 
elements  near  the  earths  surface  with  e.g.  spherical  harmonics  at  higher 
altitudes.  The  elements  are  chosen  small  near  the  earths  surface.  They 
get  larger  and  larger  with  altitude  in  agreement  with  the  potentials 
attenuation.  At  satellite  altitude  there  is  only  one  element,  or  a  small 
number  of  them. 


Theory  and  application  of  spline  functions  are  very  diversified.  There 
is  an  overlap  with  least  squares  collocation  ana  a  border  line  with 
finite  elements.  Spline  functions  were  invented  by  I.  J.  Schoenberg  and 
first  described  in  his  famous  paper  Schoenberg  Cl  946).  Since  then  they 
have  evolved  into  a  very  popular  tool  of  applied  mathematicians  as  well 
as  into  an  object  of  interest  to  theoreticians,  who  implanted  them  into 
Hilbert  spaces.  Textbooks  have  been  published,  as  for  example  Ah  I  berg 
et.  at.  C 1 967),  Boehmer  CI974). 

The  practically  minded  person  associates  with  splines  a  special 
subset  of  them,  namely  polynomial  splines.  There  is  a  widespread 
preference  for  cubic  splines.  Polynomial  splines  perform  well  in 
interpolation  problems  due  to  their  simplicity,  computational  efficiency 
smoothness  and  locality.  As  already  pointed  out  by  Schoenberg  (1946), 
one  can  construct  basis  functions  having  a  local  support. 

The  use  of  spline  functions  for  problems  of  physical  geodesy  was 
suggested  by  Davis  and  Kontis  (1970).  Me i ss I  (1 97  lb)  proposed  their  use 
for  the  representat i on  of  pointwise  known  functions  during  the  evalu¬ 
ation  of  the  explicit  integral  formulas  of  physical  geodesy,  i.e.  the 
formulas  by  Stokes,  Vening  Meinesz  and  their  refinements  due  to  Molo- 
densky.  This  proposal  was  worked  out  by  Suenke  I  (1977)  and  N'oe  (1980). 

From  the  point  of  view  of  Hilbert  space  theory  spline  functions 
are  optimal  interpolators  (or  approximators)  of  functionals.  Optimality 
relies  on  two  complementary  criteria,  namely  the  minimum  norm  property 
and  the  best  approximation  property.  Both  criteria  are  based  on  the 
choice  of  a  semi  norm.  This  choice  i s  up  to  the  user.  In  contrast  to 
general  least  squares  collocation,  certain  natural  curvature-sem i norms 
strongly  suggest  themselves  as  candidates. 

Theory  and  geodetic  use  of  splines  are  discussed  in  detail  in 
Moritz  (1978),  Lelgemann  (1980)  and  in  a  forthcoming  paper  by 
Freeden  (1981).  Collocation  with  the  use  of  parameters  can  be  inbedded 


17  - 


into  the  Hilbert  space  theory  of  splines  as  outlined  in  Boenmer  (19/4), 
chapter  4. 

Ue  shall  briefly  stress  the  point  of  v i eu  of  computational 
efficiency.  If  one  uses  generalized  splines,  as  proposed  by 
Le I gemann  (1980),  Freeden  (1981),  one  deals  with  functions  lacking  a 
local  support.  Hence  the  normal  equations  are  fu!1  This  limits  the  size 
of  the  systems  to  a  few  thousand  unknowns.  On  the  other  hand  cubic 
splines  can  be  generalized  to  2  and  3  dimensions  by  means  of  tensor 
products.  Here  basis  functions  with  local  support  are  available  Thus 
the  spline  method  competes  with  the  finite  element  method  in 
computational  efficiency.  Therefore  we  shall  discuss  this  particular 
point  in  some  deta i I 

Imagine,  for  simplicity,  a  rectangular  region  in  3-dimensional 
space  subdivided  into  box-type  elements  of  equal  size  and  shape.  The  set 
of  nodes  shall  be  identical  with  the  set  of  corners  of  the  element.  Our 
intended  use  of  the  finite  element  method  relies  on  interpolating 
functions  called  Hermite  tri-cubics.  In  any  of  the  various  boxes,  the 
function  to  be  interpolated  is  represented  by  a  polynomial  which  is  a 
cubic  in  any  one  of  the  3  variables,  provided  that  the  other  two 
variables  are  fixed.  (Considered  as  a  polynomial  in  3  variables  the 
interpolating  function  is  of  degree  9).  We  associate  with  any  node  8 
parameters  representing 
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at  this  node.  These  are  the  derivatives  of  the  potential  of  ’biaegree" 
less  than  or  equal  to  [.  By  letting  ail  parameters  having  the  value  zero 
except  for  one,  we  obtain  as  interpolating  functions  a  basis  function 
associated  with  this  particular  node  and  this  particular  parameter .  This 
basis  function  is  called  shape  function.  This  will  be  discussed  in 
detail  in  section  3.2.  Here  we  only  emphasize  that  we  have  8  basis 
functions  per  node.  Any  basis  function  is  C1  continuous  and  has  a  local 
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support  ,  "vied  to  3  adjacent  elements 

I*  „e  a  I  i  erm.o  i  :  ve : represent  our  *  :e  a  by  means  of  basis 
r one t  ens  iu:  t  up  o t  t r  -cub i c  sp !  i nes .  ue  ob tain  on i y  I  bos i s  f unc t i on 
*  or  eac.n  noae  ■  s  tne  tensor  product  3(x,y,z )  =  zhx)8(  y)3(z)  of 

one  -  a  mersional  3_solnes  3(x)  specified  ■ n  Schoenberg  ( I 94o ) ,  p.  71, 

2_ 

an  c  ca ;  ea  H (  x  )  there  Irvs  basis  function  3(x,y,z)  is  even  C 
continuous  Tnis  mexes  splines  attractive  as  compared  to  Hermite 
tr i -c^bics  however ,  it  turns  out  that  the  support  of  cny  basis  function 
covers  64  adjacent  elements 

Let  us  summarize  and  conclude  the  comparison  of  finite  elements 
and  splines  by  the  following  3  statements 

(1)  Tricubic  splines  have  basis  functions  involving  less  parameters  than 
Hermite  tricubics.  This  may  be  viewed  beneficial  On  the  other  hand  it 
also  means  that  we  have  i ess  flexibility  unless  we  decrease  the  size  of 
the  elements 

(2)  Splines  are  smoother.  This  makes  them  more  useful  in  interpolation 
problems.  However,  in  problems  of  representing  a  field  governed  by  a 
differential  equation,  we  have  an  additional  enforcer  of  smoothness 
This  was  called  the  field  contribution  to  the  normals  in  section  2.4. 

For  this  reason  splines  ore  preferred  in  pure  interpolation  problems, 
whereas  Hermite  polynomials  are  preferred  in  the  finite  element  solution 
of  field  equations  Confer  the  discussion  in  Strang-Fix  (1973),  p.  61. 

(3)  Due  to  the  larger  support  of  sp I  mes,  the  I  inear  system  ui I  I  be  less 
sparse  In  an  oversimplified  way,  we  may  talk  of  a  larger  bandwidth  as 
compared  to  a  system  resulting  from  the  use  of  Hermite  tricubics.  It 
appears  that,  whatever  may  be  gamed  by  a  decreased  number  of  parameters 
and  by  greater  smoothness  in  case  of  splines,  is  paid  for  by  a  larger 
bandwidth  slowing  down  the  elimination  procedure 
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2 . 6  Aoor ox  ? Tig te  explicite  Green’s  functions 

Frequently  the  oldest  methods  are  the  best  Hardly  ever  they  are  the 
worst  methods  computationally  Computers  were  unavailable  at  earlier 
times.  Take  Stokes/  Formula  It  yields  the  geo i da  I  undulation  at  one 
point  in  terms  of  gravity  anomalies  all  over  the  globe  It  is  hardly 
necessary  to  point  out  the  approximations  under  I y i ng  Stokes  formulc  as 
well  as  the  corrections  which  partly  make  good  For  them  The  usefulness 
of  Stokes  formula  as  well  as  of  Vening  Metnesz  formula  and  their 
refinements  is  undoubted.  Applications  are,  however,  restricted  to  areas 
of  moderately  varying  topography. 

Stokes'  and  Vening  Meinesz'  kernel  ore  explicitely  known  Green's 
functions  of  boundary  value  problems  for  the  sphere.  We  are  in  a  similar 
situation  as,  when  dealing  with  a  large  system  of  linear  equations,  an 
a-priori  known  inverse  of  the  coefficient  matrix  is  available  If  the 
boundary  value  problem  is  discretized  in  agreement  with  a  discretization 
of  the  gravity  anomalies  in  terms  of  N  block  averages,  we  obtain  indeed 
such  a  system  of  N  I  inear  equations  in  N  unknowns.  If  the  inverse  is 
known,  calculation  of  the  solution  requires  N  steps  for  one  particular 
unknown  and  N*  steps  for  a  I  I  of  them.  This  is  not  impressive  in  itself 
because  we  know  that  in  case  of  a  sparse  system  of  the  type  mentioned  in 
subsection  2.4  we  can  do  better,  namely  solve  for  all  unknowns  in  0 C N 
steps.  In  case  of  evaluating  the  discretized  Stokes  formula,  an 
important  additional  bonus  is  available,  namely  the  remote  zone  effect. 
It  implies  that  of  the  N  steps  necessary  to  evaluate  one  specific 
unknown,  many  can  be  lumped  into  comparatively  few  new  steps,  and  many 
may  even  be  omitted  altogether.  The  number  of  new  steps  to  be  carried 
out  is  a  fraction  oc.N  of  N,  where  cv.  is  viewed  as  a  fixed  constant.  The 
constancy  of  <X  is  based  on  the  following  argument.  Suppose  that  a 
certain  block  design  is  used  for  the  approximate  evaluation  of  Stokes 
integral.  Near  the  point  of  evaluation  we  use  averages  of  gravity 
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anomalies  over  the  smallest  blocks  available,  say  1 0  x  1  °  blocks  This 
corresponds  to  N  =  360*188  =  62888  At  a  moderate  distance  we  may  lump 
blocks  to  2°x2°  and  so  on  Very  distant  blocks  may  be  omitted,  in 
par  t  i  cu  I  ar  in  case  of  Vemng  Meinesz'  formula  It  we  quadruple  N, 
proceeding  from  l°xl°  blocks  to  30'x38'  blocks,  any  of  the  lumped  blocks 
in  the  above  design  is  so i : t  into  4  new  blocks  Hence  the  number  of 
steps  also  quadruples 

It  appears  that  the  effort  needed  to  calculate  the  geo i d  at  one 

2, 

point  (block  center)  is  ocN,  and  ocN  for  N  points.  Hence  the  method  is 
still  0(NZ) ,  but  the  constant  hidden  under  the  "0“ -symbol  is  very  smal I 
Finite  elements  are  0(NT),  and  consequently  asymptotically  better. 
However,  the  constant  hidden  in  OCN11)  is  large  The  break  even  point  is 
not  exactly  known  now. 

If  geo i d  or  deflection  of  the  vertical  are  needed  only  at  one 
point  or  at  a  very  small  number  of  points,  the  explicit  inverse  method 
is  the  best,  namely  0(N)  w i th  a  very  small  hidden  constant  However  it 
must  be  stressed  that  the  explicit  inverse  method  relies  on  a  special 
type  of  homogeneously  distributed  and  nonredundant  data.  There  is  no  way 
to  vary  the  weights  individually  for  the  blocks.  Additional  data  are 
difficult  to  incorporate  in  a  theoretically  satisfactory  way. 


2.7.  Exploiting  rotational  or  translational  symmetries. 

2.7  I .  Invar i ance  of  normal  equations  under  a  group  of  transformations. 


We  are  not  referring  to  a  method  that  stands  for  its  oun  as  those 
described  thus  far.  We  are  dealing  with  a  technique  that  can  be  applied 
in  conjunction  with  any  of  the  methods  described  in  subsections  2.2  to 
2.6,  provided  that  the  distribution  of  data  satisfies  certain 
requirements  which  are  rather  stringent. 


If  the  pattern  of  measurement- I ocat i ons  and  -weights  is  invariant 
with  respect  to  the  group  of  rotations  around  an  axis,  or  with  respect 
to  the  group  of  translations  in.  I,  2,  or  3  dimensions,  then  this 
symmetry  is  reflected  by  the  system  of  normal  equations  to  be  solved.  Of 
course,  proper  care  must  be  taken,  that  the  par ameter i zat i on  of  the 
potential  conforms  with  the  symmetry.  The  normal  equation  system  is 
invariant  with  respect  to  one  or  more  transformations  generating  the 
group,  provided  that  the  unknowns  are  properly  renumbered.  Let 

Gy  *  r  (2.13) 

denote  the  original  normal  equations.  Let 

X  r  H  y  (2.14) 

be  the  transformation  taking  into  account  the  translation  or  rotation 
followed  by  a  renumbering.  The  transformation  will  be  orthogonal,  i.e., 

HrH  *  I  (2.15) 

The  new  normal  equations  are 

HrG  H  y  1  HTr  (2.16) 

They  are  identical  to  the  old  ones.  Hence 

HtGH  *  G  or  GH  1  HG  (2.17) 

Two  matrices  which  commute  share  a  common  system  of  eigenvectors.  It 
follows  that  an  invariant  subspace  of  H  is  also  an  invariant  subspace  of 
G.  Invariant  subspaces  of  H  are  usually  easy  to  identify.  H  reflects 
only  the  symmetries  of  the  problem  and  is  independent  of  other 
structural  properties.  The  knowledge  of  invariant  subspaces  allows  the 
decomposition  of  the  normals  (2.13)  into  several  independent  systems  of 
sma I  I er  d i mens i on . 
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2  7  2  Outline  for  l he_c ase  of  rotational  symmetry  around  on  axis 

Assume  that  the  system  is  invariant  with  respect  to  a  rotation  around 
one  axis  by  an  angle  of 


This  case  arises,  if  5  degree  by  5  degree  mean  gravity  anomalies  are 
taken  as  measurements  and  if  the  field  is  parameterized  by  a  surface 
I ayer  u i th  cons  tanl  density  in  N  =  2592  blocks  of  size  5°x5°.  We  then 
have  k  =  360/5  =  72.  Imagine  the  parameters  (block  densities)  grouped 
according  to  longitude.  The  groups  are  numbered  according  to  increasing 
longitude.  For  a  certain  fixed  longitude,  we  imagine  a  numbering 
according  to  decreasing  latitude.  Note  that  a  rotation  by  fb  carries  all 
blocks  of  a  certain  longitude  A  over  into  blocks  of  longitude  A  +  fb  . 
Hence  a  cyclic  renumbering  of  the  groups  of  biocks  is  necessary  in  order 
to  ensure  invariance  of  the  normal  equations  under  the  transformation  H. 
H  w  i  II  be  of  the  form 


(2  18) 


The  size  of  the  diagonal  blocks  I  is  implied  by  the  number  of  blocks 
having  the  same  longitude.  Invariant  subspaces  of  H,  which  is  a 
permutation  matrix,  are  immediately  specified.  They  are  given  by  the  k 
b  I  ock  -columns  of  the  following  unitary  matrix 


•  •  •  Z  1  r-1 


7  .  —  2. 

L'  Ik  : 


(2.19) 


-r  \ 

\° . 
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w  i  th 


Zrs 


r 


r,s  -  0, 1 , K-1 


I ...  unit  matrix  .  i :  fT 


(2  20) 


A  parameter  transformat i on 


y-2z 


(2.21) 


will  decompose  the  normal  equation  system  into  k  independent  systems  of 
size  N/k .  Their  solution  will  require  an  effort  proportional  to 


(2.22) 


Since  k  =  O(VFT) ,  this  effort  is  0(N4).  Of  course,  also  the  computations 
required  to  transform  the  system  must  be  taken  into  account.  However 
here  one  may  employ  fast  Fourier  transform  techniques. 

Transformations  like  that  one  outlined  above  have  been  used  before 
in  other  disciplines.  They  have  been  used  in  geometrical  geodesy  by 
Meissl  (1969)  in  order  to  analyse  the  strength  of  regular  tr i angu I  at i on 
chains.  Their  first  use  in  physical  geodesy  is  due  to  Colombo  (1980). 


2.7  3.  The  effort  of  the  TASC  group. 

It  should  be  mentioned  that  translational  symmetry  requires  observations 
covering  an  infinite  line  or  an  infinite  plane.  If  the  domain  is 
restricted  to  a  finite  rectangle,  boundary  effects  destroy  the  symmetry. 
Nevertheless,  one  is  able  to  salvage  most  of  the  saving  encountered  in 
the  undisturbed  case.  The  methods  are  more  involved.  They  have  been 
widely  used  in  picture  processing.  Recently  an  effort  has  been  made  by 
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TASC  (-  The  Analytic  Science  Corporation)  to  analyse  geoohysical  data 
distributed  requiariy  on  an  arbitrarily  long  line  or  within  a  small 
rectangle.  Confer  Heller,  To  it  and  Thomas  (1977),  Tait  (1979) 

The  comput a t i ona i  effort  ,n  both  cases  is  proportional  to  N  log  N 
where  N  is  the  size  of  the  I  near  system  to  be  solved  Hence  N  is  equal 
to  the  number  of  unknouns  in  a  parameter  moae I  or  equa i  to  the  number  of 
measurements  in  a  collocation  model.  The  approach  of  the  TASC  group  is 
interesting.  It  rests  on  two  techniques,  namely  (!)  the 
Fast-Four i er-Transform  (due  to  Cooiey  ona  Tukey  (1965))  and  (2)  on  the 
block  decomposition  of  block  circuiant  matrices  as  outlined  above.  In 
addition  to  these  two  ingredients,  the  authors  employ  a  number  of  tricky 
maneuvers  in  order  to  tackle  the  above  mentioned  undesirable  boundary 
effect.  5y  using  a  transition  from  an  NxN  Toeplitz  matrix  to  a  2Nx2N 
block  circuiant  matrix,  and  by  employing  a  data  windou,  they  finally 
arrive  at  a  linear  system  in  transformed  space  where  a  few  diagonals 
near  the  main  diagonal  are  strongly  dominant.  After  neglecting  the 
element  outside  this  band,  and  after  adding  a  small  multiple  of  the  unit 
matrix  to  the  coefficient  matrix,  the  banded  system  is  solved  by 
Cho i esky . 

We  are  unable  to  present  the  details  in  this  report  and  refer  the 
reader  to  the  quoted  original  articles  We  only  make  the  following  4 
remarks . 

(I)  The  solution  is  approximate,  even  if  the  method  is  applied  to 
regular  data  on  a  straight  line  segment  or  in  a  plane  rectangle.  The 
errors  come  from  two  sources,  namely  (a)  the  neglect  ion  of  elements 
outside  the  band  and  (b)  the  addition  of  the  small  multiple  of  the  unit 
matrix.  The  procedure  (b)  results  in  deemphasizing  the  weights  of 
observations  near  the  boundary  of  the  region.  The  authors  give  error 
estimates  which  are  favorable.  However,  it  is  not  clear  whether  such 
favorable  error  estimates  are  available  for  other  situations  than  those 
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considered  by  the  authors.  They  test  their  method  in  predicting  gravity 
anomalies  from  geo  id  heights  by  a  collocation  procedure  Thus  they 
deduce  a  high  frequent  output  from  a  low  frequent  input.  (Cf. 

Meissl  CI97! )  for  a  discussion  of  the  frequency  content  of  various 
quantities  related  to  the  earths  disturbing  potential).  It  should  also 
be  tested  whether  deducing  a  low  frequent  output  from  a  high  frequent 
input  can  be  done  with  errors  of  the  same  small  magnitude. 

(2)  In  the  case  of  a  two  dimensional  area,  the  method  works  strictly 
speaking  only  for  a  plane  region.  Happing  a  part  of  the  sphere 
(spheroid)  onto  the  plane  causes  distance  distortions,  as  the  authors 
point  out.  However  they  do  not  point  out  that  these  distance  distortions 
interact  with  the  cover i ance  kerne  I ,  causing  it  to  fail  to  fulfill 
Laplace's  equation  any  longer.  Harmonicity  of  the  covariance  function  is 
an  inherent  assumption  in  collocation.  The  authors  map  the  systems  of 
meridians  and  parallels  onto  a  rectangular  grid  in  the  plane.  Such  a 
mapping  has  appreciable  distance  distortions.  It  is  well  known  that 
there  are  mappings  which  perform  better  in  this  respect.  Perhaps  one  of 
them  should  be  used. 

(3)  Intuitive  insight  into  the  method  can  be  gained  by  the  following 
consideration.  Consider  the  familiar  collocation  problem  for  two  random 
vectors ; 


y  =  cyx  (2.23) 

Ue  assume  that  x  and  y  are  related  to  discrete  equidistant  points  i  =  0, 
t,  .  .  . ,  n  on  the  line.  We  assume  homogeneous  processes,  hence  Cxjr  is  a 
symmetric  Toeplitz  matrix: 


(2.24) 
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We  assume  that  is  of  the  form 


C  3  C  ♦  D 

U  vv  " 


(2.25) 


where  C  is  positive  definite  Toeplitz  ond  D  is  diagonal  and  positive 
definite.  D  is  due  to  measurement  noise.  The  main  problem  is  the 
calculation  of 


7  :  (Cm  D J  ~  1  X 


i . e .  the  so  I ut i on  of 


(2.26) 


(O  D)i  *  X 


(2.27) 


As  TASC  proposes,  we  extend  C  to  a  2n  x  2n  Toeplitz  circulant  matrix  C 
where 


cK  Ok  k  <  n- 1 

0  K  *  n 

C,*„- ki  K  6  2n'1 


(2.28) 


D  will  also  be  diagonally  extended  to  D  as  we  show  in  a  moment.  We 
extend  x  by  zeroes 


*s  VO, 


(2.29) 


and  thus  obtain  the  system 

(C  +  D )  2  =  x 


(2.30) 


This  system  is  not  equivalent  to  the  earlier  one  (2.27)/  in  the  sense 
that  the  restriction  of  z  to  the  first  n  components  is  the  solution  of 
(2.27).  The  reason  for  the  failure  of  (2.30)  is  the  following  one.  The 
system  (2.30)  predicts  y  not  only  from  the  measured  x  -L  ;  i  =  0,  I,  . 
n-l,  but  also  from  the  artificially  assumed  values  xL  =  8;  i  =  n, 

2n-l.  Hence  the  prediction  is  theoretically  wrong.  However  there  is 
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sill  I  the  matrix  D  which  can  help  us  to  make  the  prediction  nearly 
correct.  This  can  be  accomplished  by  assuming  very  large  elements  of  D 
at  the  pos i t i ons  i  =  n,  .  . ,  2n- I  This  amounts  to  the  super pos i t i on  of 
very  heavy  noise  on  the  artificially  assumed  measurements  x .  =  8,  i  =  n, 
.  .  . ,  2n- 1  .  Hence  (2.30)  will  lead  to  a  near  I y  correct  pred i ct ion. 

We  now  apply  the  discrete  Fourier  transform  toward  (2  30) 
Introducing  the  2nx2n  Fourier  matrix 


,  —  ,  cr 


the  system  (2.30)  goes  over  into 


2 IT  Li  K 


(2  31) 


(2.32) 


(2.33) 


(r* A)J •  J 


r *  fhc  f 

a  -*  FmA  F 

T  -  fhx 

T  -  pH  ~ 

I  =  F  2 


Fh  is  the  Hermit i an  adjoint  of  F.  Hence  FH  is  also  equal  to  F  1. 

Tis  now  diagonal.  This  is  the  benefit  from  applying  Four.er 
transformal i on  to  a  Toeplitz  circulant  matrix.  Cf .  also  the  discussion 
around  equation  (2.20) 

A  is  a  Toeplitz  circulant.  Hence  it  appears  doubtfull  presently 
what  we  have  gained.  However  one  can  verify  that  the  subdiagonal  bands 
of  A  decrease  as  one  goes  away  from  the  main  diagonal.  This  is  perhaps 
intuitively  clear  if  one  remembers  that  the  Fourier  coefficients  of  the 


.  4m  gggjtoHfe j 
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step  f unct i on 


S(t)* 


0 

1 


O  <■  t  <  JT 
T  6  t  6  In 


ere  given  by 


— i (i-(-Vn) 

2  :r  n  l 


(2  34) 


(2.35) 


The  tapering  effect  of  the  subdiagonal  banas  can  be  made  more  pronounced 
if  the  transition  of  small  elements  of  D  in  [0,  n-l]  to  large  ones  in 
[n,  2n-|]  is  smoothed  out  somewhat.  This  causes  "data  deemphasis"  of  X;_ 
near  the  interval  ends. 

The  next  step  is  to  neglect  the  off  diagonal  bands  of  A  up  to  a 
small  number  m  (m  =  5  to  10,  or  so).  This  truncated  banded  system 
version  of  the  system  (2.32)  is  now  solved  by  Cho I esky  in  0(n)  steps. 

The  rest  can  be  accomplished  in  0(n  log  n)  steps  if  the  Fast  Fourier 
technique  is  employed. 

We  leave  it  with  this  oversimplified  picture  of  the  algorithm 
which  could  be  extended  to  the  2~d i mens i ona I  case.  It  is  interesting  to 
note  that  the  above  mentioned  "data  deemphasis  effect"  is  also 
encountered  in  the  presentation  of  the  real  TASC  algorithm. 


(4)  Shortly  before  finishing  this  report,  an  article  by 
B i tmead-Ander son  (1980)  came  to  my  attention.  The  authors  show  how  an 
nxn  Toepl  i  tz  system  can  be  solved  in  Q(n  log^n)  steps  by  a  doub I  ing 
method  The  matrix  is  subjected  to  some  mild  r estr i ct i ons,  however 
reference  is  made  to  other  work  by  Brent  et.  al.  (1980),  in  which  an 
0(n  log  n  log  n)  algorithm  is  specified  which  achieves  a  solution 
whenever  it  exists.  It  should  be  noted  that  no  approximations  are 
involved  as  they  are  in  the  work  of  TASC.  Of  course,  the  question  arises 
again  how  large  the  constant  hidden  in  the  "0"  symbol  really  is. 


Another  way  to  utilize  symmetries  was  pointed  out  by  Rauhala.  Confer 
Rauhala  (IS80)  for  details  and  further  references  Also  Snay  (1978) 
gives  an  introduction  to  Rauhala's  "array  algebra"  It  is  an  application 
of  the  concept  of  multilinear  mappings  between  tensor  spaces.  Let  X,  Y, 

Z  denote  3  vector  spaces  of  not  necessary  equal  dimension.  Consider  the 
space  T  of  tensors- 

T  --  X  G  V©  Z  (2.36) 

An  element  of  this  space  is  represented  by  a  three-dimensional  array 
t  ij  i<  -  T  can  be  viewed  as  the  linear  span  (set  of  linear  combinations)  of 
vectorial  tensor  products  x»y®z  with  x  €  X,  y  e  Y,  zeZ.  The  tensor 
generated  by  xsysz  has  elements  t^j^  =  x^Y^Zl-  Consider  three  further 
vector  spaces  X',  Y',  2'  of  arbitrary  dimensions  I',  J',  K'.  Form  the 
tensor  product  T'  =  X'®Y'®Z'.  Let  A,  B,  C  be  linear  operators 

x’*  Ax 

y  *  8y  (2.37) 

Z':CZ 

Define  a  linear  map  T-»T'  in  the  following  way.  Let  the  image  of  x®y®z 
be  Ax®Ay®Az.  Extend  the  domain  from  the  set  of  tensor  products  to  all 
of  T  by  means  of  linearity.  Thus  a  map  A®3oC  from  T-*T'  is  obtained. 
Suppose,  temporarily,  that  any  of  the  maps  A,  B,  C  is  invertible.  It  is 
easily  proved  that 

(A©8  ®  C)'1  *  A*1®  B‘/|  ®  C’1  (2.38) 

Here  the  benefit  from  array  algebra  becomes  transparent:  instead  of 
inverting  a  huge  matrix  of  size  (IJK)#(IJK),  one  inverts  3  matrices  of 
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sizes  1*1,  J*J,  K*K . 

Assume  nou  I'  >  I,  J'  >  J,  K'  >  K  and  consider  the  least  squares 
prob I em 

t‘ ♦  V  3  (A®B®C)  t  (2.39) 


with  t  6  T,  t'e  T',  and  veT'  denoting  the  residuals.  The  3  least  squares 
prob ! ems 

x  *  v,  *  Ax 

y‘  *  v,  -  By  C2.40) 

2*  *  vt  --  Cz 

are  solved  by  the  pseudo- i nverses 

X  *  A+  y' 

y  *  BV  (2.41) 

2  --  C'z' 

If  the  rank  of  A  equals  the  number  of  its  columns,  then 

A+  3  (ArA)'1Ar  (2  42) 


and  similarly  for  B*,  C"*". 

Rauhala  shows  that  the  pseudo- i  nverse  of  A®B®C 


is  given  by 


(A  ®  B  ®  C)+  =  A  ®  B  ®  C 


(2.43) 


A  proof  follows  eas i  ly  from  geometric  reasons  based  on 
“null-space"  considerations.  We  do  not  give  a  complete 
because  it  requires  a  number  of  formal  definitions.  We 
the  following  facts : 


"range-space’  and 
proof  here, 
mere  I y  ment i on 


(a)  The  matrix  A  maps  its  domain  space  one  to  one  onto 


its  range  space. 
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The  orthocomplement  of  the  domain  space  is  the  nu I i -space 

(b)  The  pseudo i nverse  maps  the  range  space  inversi ly  back  onto  the 
domain  space.  It  maps  the  orthocomplement  of  the  range  space  onto  zero. 

(c)  Domain  space  and  range  space  of  A®B®C  are  the  tensor  products  of 
domain  and  range  spaces  of  A,  B,  C.  Hence  (2  43)  is  essential ly  reduced 
to  C2.38). 

Thus  the  least  squares  problem  (2.39)  is  solved  by 

t  ■  (A*®  B*@  O  t'  «■«« 

Let  this  outline  be  enough.  We  just  mention  that  generalization  to 
tensor  products  of  arbitrarily  many  factors  are  immediate.  Rauha I  a  will 
also  forgive  that  I  did  not  use  the  most  generalized  inverses  he  has 
ever  invented. 

Applications  of  array  algebra  are  restricted  to  gr i dded  problems 
defined  on  regions  being  of  the  box-type.  The  grids  must  be  rectangular, 
however  the  spacing  between  gridlines  may  very.  The  question  still 
remains  how  familiar  problems  of  physical  geodesy  are  transformed  into 
problems  of  array-a I gebr a .  Rauhala  states  that  this  cannot  be  done 
without  some  “cheating  a  la  Gordian  knot".  The  cheating  may  perhaps  be 
comparable  to  that  one  encountered  during  the  transformation  of  a 
problem  formulated  for  a  small  spherical  rectangle  into  a  translation 
invariant  problem  defined  over  a  plane  rectangle. 

Let  us  discuss  the  computational  effort  for  the  case  of  two  vector 
spaces  X,  Y  of  equal  dimension  n.  The  tensor  product  T  =  X®Y  has 
dimension  N  =  n1.  The  effort  to  naively  solve  an  N*N  system  requires 
0 C N 3 )  =  0(n6)  steps.  The  effort  to  invert  2  matrices  of  nxn  is  Q(n3)  = 
O(N^) .  It  may  be  shown  that  a  solution  of  the  system  utilizing  the 
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i 

decomposition  A®B  ccn  also  be  done  in  OCNO  steps.  This  is 
asymptotically  the  same  effort  as  if  a  general  sparse  system  resulting 
from  a  2~d i mens i ona I  layout  is  solved  by  the  nested  dissection  method. 
Confer  the  discussion  in  section  2.4,  in  particular  the  first  of  the  2 
remarks  given  at  the  end  of  subsection  2.4 


I? 
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g  i  ven 

f  (0)  '  C„  fW  ■  c,0 

{'(o)  •  cQ1  f'd)  ■ 


f i nd  a  cubic  polynomial  i nterpo i at i ng  these  values. 

Z  cLiY£j(*)  ;  0  ±  x  ±1 

c.j.O 

Ue  also  introduce  the  derivatives 

i  k 


the  solution  is 

(3  5) 


rkiJ  w 


Til  (*)  •  K--  o,  i,  i 


(3.6) 


Remark :  The  coefficients  of 


rKij  M  ■  E  cKiit  xe 


£.0 


KtJ  £ 


(3.7) 


are  stored  in  the  4-d i mens i ona i  ar.  ay  PSC(K,I,J,L)  during  execution  of 
our  computer  programs  described  in  chapter  6 


3  12.  Interval  of  length  h 
Let 


Yij  (x;h)  ’  h J  TcjI'k)  ;  0  £  X  i,  h  (3.8) 

These  functions  solve  the  above  interpolation  problem  with  data  given  at 
x  =  0  and  x  =  h.  Again  we  take  derivatives  of  these  basis  functions: 

■Wy;M '  its 


k  •  O.i.  2. 


(3.9) 
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1 


Obv 


i  ous I y : 


Tkijt*; h)  *  hJ'kYkij 


(3. 10) 


3  I  3.  Bicubic  polynomiol  in  a  rectangle  with  sides  a,  b 


Take 


'Y' 


LxLyJx)y 


(y.y;a,b)  *  y.  ■  (x;a)  ■  v\  .  (y;b) 

*  **  •  Ly Jy 


o  ±  x  £  a 
0  4  y  6.  b 


(3.11) 


Note:  the  first  two  indices  refer  to  the  location,  the  second  pair 
refers  to  the  degree  of  derivatives! 

The  functions  (3.1!)  solve  the  following  interpolation  problem: 

G  i  ven 


f(O.O)  =  coooo 

f(o.b)  •  C„co 

fy  (Q®)  *  C  oool 

fy  (O.b)  •  C  OIOI 

fJO.O)  •  c00„0 

*  C-oeil 

C,111 

find  a  bicubic  polynomial 

3  3 

f(y,y) 3  Z  Z 

xr>s 

r-.o  j«o 


(3. 13) 


i nterpo ! at ; ng  these  data  The  solution  is 


f  (y.y)  •• 


Z  . 

lx 


'V'. 

*  ‘■x  uy J >r  Jy 


(x'.y.-a.b) 


O  i  x  i  o, 

0  h  y  £  b 


(3.14) 


.  4 .  Tr i cub  i  c 


box  w i th 


The  extension  to  3  dimensions  is  obvious:  In  the  interpolation  problem 
we  prescribe  the  nodal  derivatives  shown  in  table  3.1  at  the  8  corners 
of  a  box  with  s i de- I ength's  a,  b,  c 


j*  i  h  i 

OOO 

ooi 

010 

01 1 

1 00  101 

1 10 

111 

Derivative 

f 

h 

fy 

fyz 

fx  fx  i 

fvy 

fx yi 

Table  3.1 

Numbering  of  nodal  derivatives  in  3  dimensions 


cont inu i t\ 


I  ement  bounder i 


Consider  an  n  =  1 ,  2  or  3-dimensional  region  divided  into  elements.  For 
n  =  1  the  elements  are  intervals,  for  n  =  2,  3  the  elements  are 
rectangular  boxes.  For  n  =  2,  3  we  assume  that  a  corner  of  an  element 
does  not  touch  the  interior  of  a  face  (boundary  segment,  boundary 
rectangle)  of  an  adjacent  element.  Otherwise,  the  elments  need  not  be 
of  the  same  size.  Take  as  an  example  the  partition  of  a  region  in  R1 
shown  in  figure  3.2. 
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F i gure  3 . 2 

Sample  element  partition  in  2  dimensions 

Assume  that  the  Hermite  nodal  values  are  prescribed  at  the  corners.  Then 
a  function  may  be  interpolated  into  any  of  the  elements.  It  is  an 
n-cubic  polynomial  there.  This  function  is  continuous  and  has 
continuous  first  derivatives  everywhere.  Such  a  function  is  said  to 
belong  to  the  class  C*. 

The  proof  is  easy  but  not  entirely  obvious.  Let  us  sketch  it  for 
n  =  2.  Take  any  segment  separating  two  elements,  e  g.  the  line 

segment  A-B  in  figure  3.2.  The  limits  of  f(x,y)  from  the  top  and  bottom 
elements  are  two  cubic  polynomials  f^pCx),  fsorrowCx)  in  x.  However  such 
polynomials  are  completely  determined  by  the  values  f(A),  f*(A),  f(B), 
fjfCB)  which  are  common  to  both  polynomials.  Hence  the  polynomials  must 
coincide.  The  reader  may  wish  to  extend  the  argument  to  the  continuity 
of  fy  across  the  segment  A-B.  It  becomes  transparent  why  the  mixed 
derivatives  f*y  are  needed  at  the  nodes. 
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3. I .6.  Compound  elements. 

We  will  encounter  functions  which  are  rough  in  some  areas  and  smooth  in 
others.  In  order  to  rep-esent  them  properly  and  economically  we  like  to 
be  able  to  change  the  size  of  the  elements  in  a  way  that  is  more 
flexible  than  that  indicated  in  figure  3.2.  In  order  to  achieve  this,  we 
must  sacr i f i ce  something,  namely  the  simplicity  of  the  elements.  It  will 
be  sufficient  to  outline  the  procedure  in  R1  and  to  consider  the 
"compound  element"  shown  in  figure  3.3 


Figure  3.3.  Sample  of  a  compound  element 

The  compound  element  shown  has  n  =  7  nodes.  At  each  node  i  four  nodal 
parameters  are  prescribed.  They  are  fCi),  fy(i),  f*(i),  f C i 3 .  If  a 
point  u  is  to  be  interpolated  which  is  situated  in  the  upper  quad,  we 
just  take  formula  (3.13)  specified  above.  Note  that  also  the  artificial 
node  h  may  be  viewed  as  a  node  of  the  upper  quad.  Its  nodal  values  fCh), 
fyCh),  fr(h),  f xy C h 7  are  interpolated  linearly  from  those  of  i  =  4  and 
i  =  5  alone  and  do  not  depend  on  those  of  i  =6,  i  =7.  Having  nodal 
parameters  in  h,  it  poses  no  difficulties  to  interpolate  points  in  the 
lower  quads.  The  interpo I  at  ion  formula  for  any  point  x  in  the  compound 
quad  may  now  be  written  as: 

f(*)  *  t  i. 

'  t-l  j*f  J  J 


(3  15) 
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The  index  i  refers  now  to  the  nodes,  the  index  j  to  the  nodal 
derivative  The  basis  functions  s"Lj(x,y)  of  the  compound  quad  are 
piecewise  cubic  polynomials.  They  are  composed  of  cubic  polynomials 
having  as  domain  the  3  subquads  making  up  the  compound  quad.  Note  that 
the  artificial  node  h  does  not  enter  the  interpolation  formula.  Its 
nodal  parameters  are  I  inear  functions  of  those  at  the  genuine  nodes,  and 
have  thus  been  eliminated.  Formula  (3.15)  holds  also  for  simple  quads, 
in  which  case  the  cr^  (x,y)  are  just  the  YlwLyJ.jy  (x,  y)  'n  a  different 
notat i on . 

C”1  continuity  within  the  compound  quad  is  obvious.  It  is  further 
obvious  that  C1  continuity  holds  within  a  region  partitioned  into  simple 
and  compound  quads  in  a  way  that  any  node  of  a  quad  is  shared  by  a  I  I 
neighbouring  quads.  Figure  3.4  gives  an  example  of  such  a  region. 


using  simple  and  compound  elements 
Of  course,  alternative  shapes  of  compound  elements  can  be  designed. 
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Figure  3.5  shows  some  of  them. 


F i gur e  3 . 5 

Other  examples  of  compound  elements 


The  idea  is  always  the  same;  decompose  the  elements  into  subelements  of 
simple,  i.e.  rectangular  shape.  Make  a  choice  of  nodes  which  you  like  to 
retain  in  the  f  inal  compound  element.  The  nodal  parameters  at  the  other 
nodes  must  be  linearly  expressible  in  terms  of  the  parameters  of  the 
retained  nodes  by  linear  interpolation  (not  extrapolation!).  After 
specifying  nodal  parameters  at  the  retained  nodes,  interpolation  of  any 
location  is  done  by  a  formula  like  (3.15).  The  basis  funct  i  ons  <y(x,  y) 
are  obtained  by  specifying  parameter  values  equal  to  zero  except  for  one 
parameter  where  a  value  of  !  is  specified. 


3.2.  Shape  functions. 

For  presentational  purposes  it  is  useful  to  introduce  shape  functions. 
Consider  a  region  subdivided  into  finite  elements  as  for  example  that 
one  shown  in  figure  3.4.  Label  the  nodes  in  some  way  by  i  =  1,2,  .  .  .  At 
each  node,  label  the  nodal  parameters  (e.g.  f,  fy,  f*,  fXy)  by  i  = 

I, 2, 3, 4.  For  each  pair  (i,j)  consider  a  function  S  ij  (x  , y)  which  has  zero 
nodal  values  at  all  nodes  i '  4-  i .  Furthermore,  at  node  i,  all  nodal 
parameters  j'  +  j  are  also  zero.  Si,j(x,y)  will  be  nonzero  only  in 
elements  adjacent  to  the  node  i  This  feature  is  a  great  advantage 


because  precisely  the  locality  of  the  shape  functions  $ij(x,y)  is 
responsible  for  the  sparsity  of  the  normal  equations  to  be  der i ved 
later.  In  a  particular  element  adjacent  to  node  i,  Sij  (x,y)  of  course 
coincides  with  the  local  functions  <T;j(x, y)  introduced  earlier  (cf.  equ. 
(3.15)).  Hence  the  shape  functions  S;.j(x,y)  are  nothing  new.  In  the 
later  developments  we  will  mainly  deal  with  their  fragments,  the  element 
related  &ij  (x, y) ' s .  However,  many  interrelationships  are  more  clearly 
explained,  if  the  globally  defined  functions  Sij(x,y)  are  used. 

3  3.  Contribution  of  the  field  to  the  normal  equations. 

In  conventional  least  squares  setups,  the  normal  equations  are  formed 
from  observations.  Observations  will  also  contribute  to  the  normal 
equations  in  our  case.  However  there  will  be  another  contribution.  The 
harmonicity  of  the  field  is  not  automatically  implied  by  the  Hermite 
cubic  representation.  Complete  harmonicity  is  practically  incompatible 
with  this  field  representat i on .  All  that  can  be  done  is  an  approximate 
fulfillment  of  Laplace's  equation.  This  will  be  achieved  by  minimizing 
the  integral  over  the  squared  Laplacean  of  the  field  This  integral  will 
give  the  additional  contribution  to  the  normals  as  announced  ear  I i er . 


3.3.1.  Reasons  for  excluding  the  Ritz  method. 

Our  least  squares  approach  may  be  called  an  old  fashioned  one.  At  least 
this  is  indicated  in  current  treatises  on  finite  elements  (cf.  Strang; 

Fix  (1973),  pp .  133  to  134).  In  these  treatises,  attention  is  focused  on 
more  modern  principles  due  to  Ritz  and  Trefftz  and  generalizations  of 
them  published  in  the  last  few  decades.  Confer  Oden  (1979).  Let  us 
explain  why  we  do  not  propose  the  Ritz  principle.  (I  tried  it  in 
numerical  experiments,  but  it  gave  poor  results!).  The  Ritz  principle 

i 
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successfully  used  in  the  following  typical  S'tuation 
Find  a  solution  of 


AV  *  0  in  8 

V  1  g  on  <98 


(3.16) 


This  is  the  familiar  Dirichlet  problem.  However  not  3 i r i ch I et i c i ty  is 
the  point.  We  could  have  used  Neumanns  boundary  conditions  as  well,  or 
even  a  mixture  of  both.  The  decisive  point  is  that  fixed  boundary  values 
are  prescribed.  The  Ritz  principle  replaces  the  above  proolem  by  a 
var i at i ona I  one : 

Find  a  solution  of 


\  f  I  9ra°L  V  |Z  dB  »  Min. 


(3.175 


subject  to 


V  9  g  on  dB 


(3. 18) 


The  variational  formula  is  slightly  more  complicated  in  case  of  Neumanns 
boundary  condition,  but  this  is  irrelevant  presently. 

The  next  step  is  to  replace  V  in  (3.17)  by  its  finite  element 
representation,  i.  e.  by 


V  ■  £  V„.(*.y) 


(3.19) 


Here  Vy  are  the  unknown  nodal  parameters  and  Sij(x,y)  are  the  known 
shape  functions.  The  functional  (integral)  to  be  minimized  becomes  a 
quadratic  function  in  the  unknowns  V^.  The  fulfillment  of  the  boundary 
conditions  V  =  g  at  98  can  not  be  postulated  in  a  strict  sense.  One  has 
to  be  satisfied  that  V  =  g  at  certain  points  of  the  boundary 
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and  perhaps  also  that  some  derivatives  of  V  and  g  along  the  boundary 
coincide  at  these  points  The  boundary  conditions  thus  yield  a  linear 
set  of  constraints  for  the  unknowns  V  LJ .  Minimization  of  a  quadratic 
functional  subject  to  linear  side  constraints  is  a  standard  problem 
which  leads  to  the  familiar  linear  normal  equations  whose  solution  are 
the,  V  ij  . 

Consider 'a  sequence  of  partitions  of  B  into  finite  elements  such 
that  the  diameter  of  the  largest  element  goes  to  zero.  It  is  also 
required  that  the  shape  of  the  elements  is  not  too  badly  distorted  as 
the  diameter  tends  to  zero.  In  treatises  on  finite  elements  it  is  proved 
that  under  fairly  general  conditions  the  finite  elment  solution 
converges  to  the  exact  solution  of  the  original  problem.  The  main 
advantage  of  the  Ritz  method  over  the  more  primitive  method  indicated 
above  and  to  be  described  in  detail  below,  namely  the  method  of 
minimizing 


(3  20) 


subject  to 


V  --  5  at  dB  (3.21 ) 

is  the  following  one:  The  Ritz  method  involves  a  functional  defined  in 
terms  of  first  derivatives.  As  a  consequence,  one  may  use  shape 
functions  Sij(x,y)  which  are  simpler  than  those  required  for  the  other 
method.  The  functional  (3.20)  is  defined  in  terms  of  second  derivatives 
It  is  shown  in  the  literature,  that  our  piecewise  cubic  polynomials, 
which  are  C1  continuous  across  element  boundaries,  are  an  admissible  set 
of  trial  functions  for  the  least  squares  problem  (3.20,  3.21).  However 
the  Ritz  problem  can  be  treated  successfully  with  trial  functions  which 
are  only  piecewise  bilinear  and  just  C  continuous  across  element 
boundaries.  Such  functions  require  only  one  nodal  parameter,  namely  the 


function  value  f  at.  the  node  This  holds  in  1,  2  on  3  dimensions.  Recall 
that  our  C  continuous  shape  function  require  2“*  nodal  parameters,  i.e. 
2,  4  or  8  depending  on  the  dimension  d.  The  decrease  in  number  of 
parameters  offers  one  advantage  although  it  must  be  ba I anced  against  the 
need  to  use  smaller  elements  because  of  the  more  primitive  nature  of  the 
shape  functions.  In  any  case,  the  computer  programs  turn  out  to  be  much 
simpler  if  bilinear  shape  functions  are  used. 

But  why  does  the  Ritz  method  not  work  in  our  geodetic  environment? 
In  geodesy  the  boundary  values  are  the  result  of  measurements. 
Measurements  are  not  performed  everywhere  and  they  are  subject  to 
observation  errors.  In  some  areas  the  measurements  are  redundant,  e.g. 
there  may  be  measurements  of  V  as  well  as  of  some  components  of  the 
gradient  of  V  In  other  areas  the  measurements  may  be  sparse.  There  may 
also  be  measurements  in  the  interior  of  B.  The  constraints  in  our 
minimization  problem  become  weighted  constraints,  so  to  speak.  A 
problem  of  balancing  the  weights  arises.  The  normal  equations  are  the 
sum  of  two  contributions,  namely  that  one  from  the  minimization  of  the 
functional,  and  that  one  from  the  observation  equations.  Symbolically 

(PlG1  *  P2Gz)  x  *  r  (3  22) 

Suppose  that  the  functional  to  be  minimized  is  the  R i tz-funct i ona I ,  i.e. 

1  f  |  9rad  V  pets  (3.23) 

3 

If  we  choose  p,  large  and  p4  small,  then  the  gradient  of  V  will  be  made 
small  at  the  cost  of  large  residuals  at  the  observat i ons .  In  the  limit 
p  we  get  a  constant  V.  If  we  choose  px  large  in  comparison  to  p., 

we  treat  the  observation  equations  practically  as  constraints.  The 
minimization  procedure  will  try  to  match  the  observations  exactly.  This 
may  lead  to  absurd  results  in  case  of  redundant  observations.  Even  in 


-  45  - 


case  of  nonredundancy,  the  resulting  function  V  may  be  far  off  from  thct 
one  minimizing  the  gradient.  Hence  harmonicity  is  not  ensured.  Another 
disadvantage  is  encountered  in  subsequent  error  propagation  studies  If 
one  tries  to  propagate  observation  errors  to  some  desired  quantity,  eg. 
a  gradient  at  satellite  altitude,  the  resulting  error  will  reflect  the 
small  p.,  rather  than  the  large  pa.  This  is  extremely  undesirable, 
because  the  propagated  error  is  then  mainly  due  to  a  "poorly  observed" 
gradient  rather  than  to  the  observation  errors.  These  disadvantages  are 
to  a  large  extent  avoided  if  we  choose  the  minimizing  functional  as 


(3.24) 


From  the  above  discussion  it  is  clear  that  weights  p1  should  be  rather 
large  in  comparison  to  pt.  p.,  weighs  now  the  failure  of  V  to  be 
harmonic,  whereas  in  the  earlier  case  it  weighted  the  failure  of  V  to  be 
a  constant  function.  There  is  a  lot  of  harmonic  functions  which  deviate 
considerably  from  a  constant  function. 

If  we  minimize  the  functional  (3.24),  there  will  be  deviations 
of  AV  from  the  zero  function,  but  these  deviations  will  be  small.  Hence 
the  weights  p1  can  be  chosen' large  in  comparison  to  pr.  Errors 
propagated  after  adjustment  will  therefore  mainly  reflect  the  errors  due 
to  the  smaller  weights  pt  which  belong  to  lh^  observations. 


3.3-2.  The  least  squares  contribution  of  the  field. 

After  discussing  the  reasons  for  choosing  the  functional  (3.24)  for  the 
purpose  to  ensure  approximate  harmonicity,  we  turn  to  the  technicalities 
of  calculating  the  contribution  of  this  functional  to  the  normal 
equations.  Let  us  first  outline  the  procedure  on  hand  of  the  shape 
functions  Sij(p).  As  explained  earlier,  the  index  i  refers  to  a  node 
while  the  index  j  refers  to  one  of  the  parameters  at  this  node  The 
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argument  p  may  now  be  viewed  as  !,  2  or  3  dimensional  J.n  case  ot 
detailed  discussions  we  shali  mainly  stress  the  case  of  2~d i mens i ona I 
poiar  coordinates,  where  we  write  SLJ(r,y?)  instead  of  SLj(p) 

The  potential  is  represented  cs 

V  ‘  ^  K]  (P)  (3  25) 

tj 

Its  Lap  I acean  is 

AV  1  Z  V--  A5L-(p)  0  26) 

Lj  J 

If  2-dimensional  polar  coordinates  are  used,  the  Lap  I acean  is 

&V  -*  Vrr+jVr  +  Ji  (3  27) 

In  order  to  have  something  specific  in  mind,  the  reader  may  imagine  the 
region  3  in  the  form  of  a  circular  ring  In  polar  coordinates,  the  ring 
becomes  a  rectangular  region.  It  may  be  subdivided  into  simole  and 
compound  elements  as  shown  in  figure  3.6(a)  The  subdivision  shown  there 
corresponds  to  a  subdivision  of  the  ring  which  is  shown  in  figure 
3.6(b). 


Figure  3.6. (a).  Partition  of  a  circular  ring 
represented  in  poiar  coordinates 


Figure  3.6. Cb) .  Partition  of  a  circular  ring 
represented  in  cartesian  coordinates 

Such  a  subdivision  may  be  chosen  if  the  field  is  anticipated  to  be  more 
detailed  in  the  vicinity  of  the  inner  boundary.  Substituting  the 
Laplacean  (3.27)  into  the  functional  (3.24)  gives 


T  f  (  ZV^AStj.fp))1  dB(p) 

Jv>  LJ 


(3.28) 


«*  »Vk « w 
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Minim, zing  with  respect  to  the  Vtj  leads  to  the  normal  equations 


9 g 

3V:] 


--  0 


i  . .  nodes 
j  .  .  par ame tens/ node 


or 


£ 

a* 


f  ( p)  &SLy  (p)  oL8(p) 
& 


v,y 


i  ...  nodes 
j  ...  parameters/node 


(3.29) 


(3.30) 


Up  to  a  multiplication  of  these  equations  by  the  weight  p1t  we  get  the 
contribution  of  the  functional  (3.24)  to  the  normal  equations,  ^ote  that 
for  a  particular  equation  labeled  (i,j),  there  are  only  a  few  nonzero 
coefficients.  If  we  write  the  equations  as 


then 


»  0 


(3  31) 


^  ij;  I'j' 


0 


(3.32) 


unless  the  nodes  i  and  i'  are  "neighbored".  Two  nodes  are  neighbored  if 
there  is  an  element  in  which  both  nodes  participate.  The  sparse 
structure  of  the  normal  equations  becomes  visible  now. 


3.4.  Detailed  instructions  for  computer  implementation 

If  one  works  with  a  computer,  the  above  outl i ned  procedure  is  not 
recommended.  It  is  preferable  to  decompose  the  total  "energy*  E 


into  contributions  of  the  individual  elements  now  denoted  BK. 


E  ■  E  Ek 

K  K 


(3  33) 


w  i  th 


T f  [Ev/,.,a^.,(p;]V8(p) 

B  L.  *" 


(3.34) 


In  agreement  with  the  earlier  notation  trLj  (p)  for  the  e  I  ement- i  nterna  I 
basis  functions,  this  may  be  written  as 


l/.j,  Aflr£iyi  (  p) 


dL3(P) 


(3.35) 


The  summation  needs  to  be  extended  only  over  nodes  participating  in  the 
element.  Ue  may  form  the  'partial  normals"  due  to  the  contribution  of 
element  k: 


(3.36) 


or 


lJ 


M  I /  *  0 


(3.37) 


If  there  are  I  nodes  participating  in  the  element,  and  if  J  is  the 
number  of  parameters  per  node,  the  above  is  a  symmetric  system  of  IJ 
equat i ons . 

The  normal  equations  (3.31)  for  the  entire  region  are  obtained  by 
summing  all  partial  normal  equations.  The  situation  is  similar  to  a 
network  adjustment  where  partial  normals  may  be  formed  for  subsets 
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of  the  measurements;  a  subset  may  e.g.  be  given  by  a  round  of  directions 
taken  at  one  station.  The  par 1 1 a  I  normals  are  then  combined  to  yield  the 
normals  of  the  entire  network. 


3.4.1.  The  partial  normals  of  a  simple  auca. 


We  return  to  the  notation  of  section  3.1.2.  for  the  basis  functions  in  a 
simple  quad  writing: 


V(ry)  --  z  V  ,  n  ■; '  (r-n, ;  r,-  re)-  r  ,  {</-% „■  %-?0) 


(3.38) 


The  simple  quad  is  assumed  to  extend  over  the  region 

6  r  4  n, 

&  *  y  & 


(3.39) 


Our  notation  appears  complicated  now.  However  we  are  trying  to  give 
detailed  instruction  for  an  efficient  computer  implementation  rather 
than  trying  to  please  the  casual  reader  with  some  slick  and  polished 
notation  which  suppresses  the  nasty  details. 

We  form  the  Lapiacean  (confer  (3.9)  on  the  notation  for  the 
derivatives  of  the  y^s) 


Avln ’  cavVv  'Ww  InwA-*: r.w 

*  T  Viir'Jr'  (r'r''' r>  'r‘>  r.  L,y  ( f-% ;  y.-%)  t 

*  7?  roWr.  lr-r.;r.-r,,) (y-y.; 


(3.40) 


r 
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Let  us  abbreviate  this  for  the  moment  as 


AV  (r.cf)  =  Y. 


vww  Sfl"l’(r,hw(^ 


(3.41) 


In  order  to  form  the  normals  (3.36,  3.37)  ue  have  to  evaluate  all 
integral s 


I 


KK' 


•t  t-y  JrJf 


W,;W*V  ■  / /fZ(rlf^(r,Kh(^KiM  rdrd* 


r.r0  yy. 


(3.42) 


The  coefficients  of  the  partial  normals 


,•  i^r  •  i"  i  ^  *r  Hr  J  r jtr !  Lr'  i»'irV  'Jr'Jy1 


irVJr'jV  1r  HpJrJif  -*  «r'  Lf  Jr'j-/' 

are  then  given  by 


(3.43) 


y  -j-kic' 

^  i-r^iriy  !  ir’Vjr'V  K,k»1  ^r' 


y'jr'Jy1 


(3.44) 


It  is  a  great  computational  advantage  that  the  above  integrals  decompose 
into  one-dimensional  integrals 


-j-  KK' 

£rhfVi<f  y  y  Jr  Jf 


—  KK'  -y-  K*c' 

LrirWr’  yiy  W 


(3.45) 


w  i  th 


r, 

r  k*'  --  f  f Ck)  (r)  •  f  ,  W  ^ 

“*■  irjr  W  J  ll-rJr  Lr  Jr 


(3  46) 


r.r. 


* 


(3.47) 


yy 


It  is  suf f : c  lent  to  outline  the  further  procedure  on  hand  of  one  of 
these  integrals.  Let  us  take  as  an  example 


cirW  'friirir ir-r.;r,-r0)  T^jAr-r^r.-r.)^ 


(3.485 


We  subst i tute 


r-  r^ytr.-rj 


(3.49) 


i 

'/wjtr-rjsvrj  r„vv(?Cr--r”);':'i:)  r.r/(l!-r.)  ^  0.50) 


By  (3.18)  in  section  3.1.2.  this  further  equal: 


The  funct  i  ons  Yt,irjr^?^  are  now  derivatives  of  the  basis  functions  for  the 
unit  interval  Confer  section  31.1.,  equation  (3.6)  They  are  cubic 
polynomials  whose  coefficients  are  stored  in  the  array  PSC(K,I,J,L) 
during  the  computer  runs. 

According  to  the  outlined  rule,  we  get  for  any  of  the  integrals 


J.  kk' 

trjrlr’ir' 


an  expression  of  the  form 


(3.52) 


—KK‘  .  .  Jr  +Jr  *  dr(K')+  1 

IiWrW' 


■ jTi4*Kr\r (?)  TirOrtirir'  (P>  (fi  .f(r,-r0!)‘<K,'e<K,,',d5j 


(3  53) 
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ana  I ogous i y 


LWhWh' 


■fr^w  iy  jv ^  r  w  m 

o 


(3.54) 


Thereby 


CJ  3 


y  -  fo 
if*-  *© 


(3  55) 


The  following  table  lists  the  degrees  dr(k),  dy(k)  and  the  exponents 
e(k) 


k  dP  dy  e 


12  0  0 

2  1  0-1 

3  0  2  -2 

Table  3.2 

Listing  dr(k),  dy(k),  e(k);  k=l,2,3 


All  these  integrals  are  extended  over  weighted  products  of  the  basis 
functions  and  their  derivatives  for  the  unit  interval.  These  basis 
functions  are  cubic  polynomials  in  one  variable.  The  first  and  second 
derivatives  are  quadratic  and  linear  polynomials  respectively. 

Hence  all  the  integrals  can  be  assembled  with  the  help  of  the 
PSC-table  and  the  following  integrals  over  weighted  monomials: 

A 


M 


'JxHr0<-  (n,-r„)x)taix 


O 


0  i  $  $  6 

-3  i  t  $  1 


(3.56) 
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Summarizing,  we  obtain  the  following  procedure  to  evaluate 
Scr4^jrjy;  ir'  iy'jf'jy1 

for  a  simple  element  fo  £  f  i  ;  i/o  -  !f  s  if- f 

[ I ]  Evaluate  the  weighted  integrals  over  monomials  according  to  (3.56) 
and  tabu  late  them  for  0  4  s  4  6,  -3  4  t  4  I 


[2j  Using  the  PSC-tables,  whose  coeff  i  ec  i  ents  are  denoted  c^rj(.ci. 
equation  (3.7)1),  and  the  table  3.2  for  d(k),  e(k),  evaluate  integrals 
over  products  of  normalized  basis  functions 

[ TjfM  UV  <?>  t  r.  *  (  r,  -  r0)  y]  dp  - 


3  3  KK; 

§  ^  Cd-r(Ky''rJct0J.r(«>Ur,jr'r^t*e'l  eOti+eW  +  1  *  ^rirW 


J  ri,(KiiyJv(o)r<t,(xvvjv(^  ^ 


(3.57) 


3  3  KK; 

^  C<VK> ufjyC<A^KV Uf'jy'e'^t+e', o  Wh'-'iW 

[3]  According  to  equations  (3,53),  (3.54)  evaluate; 


j r  ♦  J r *  -  ir(K>  -  -,KK’ 

‘rJr'-rjr 


T  .  .  *  f  r  -  r-  V  wr  r'"’  .  'p'' 

trlrtrUV-  ^  1  f 


T**’  ,  ,  \  jvtjj-  al^CK)-  c4(tC)  M  —  KK1 


(3.58) 


0  i  i  i  L  ‘  i  1  L  i  L  ’  i  1  6  1 
rJr  r  Jr  .  ‘•yJy  y  Jy  ' 


and  tabulate  these  integrals 
[41  Assemble 


«P  i  / 

X"  x KK  x KK 

k  K  ;  1  ‘•rjr  <VjV  w ■ if  tt‘jV ' 


(3.59) 


Thus  the  partial  normals  for  a  simple  element  are  evaluated! 

Remark •  There  are  some  further  computational  shortcuts.  Because  of 
symmetries,  some  of  the 


—  kjt1 


j-  KK1 

“■mWjV 


are  the  same.  Such  symmetries  could  be  exploited,  but  we  have  refrained 
from  doing  so  in  our  numerical  experiments. 

Remark :  The  mu  1 1 i - i ndex-notat i on  is  convenient  for  calculating  the 
coefficients  9irivjrjv-  ir'iT‘jr'jv'  ■  ^or  subsequent  calculations  it  is 
preferable  to  switch  to  a  two-dimensional  array  gai.  The  mu  I  t  i  -  i  ndexes  ir 
ty imply  a  node  numbering  n  =  1,2, 3, 4  as  shown  in  figure  3.7  and  given 
by 

n  =  2iP  +  iy  +  1  (3.60) 


Figure  3.7.  Numbering  of  nodes  in  a  simple  quad 


i 
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The  parameter  numbering  per  node  has  been  explained  above  repeatedly,  It 
is  given  by 

q  =  2jr  +  jy,  +  I  <3.61) 

and  shown  in  table  3.3 

q  jr  j  Math,  symbol 

10  0  V 

2  0!  Vy 

3  10  VP 

411  Vry 

Table  3.3.  Parameter  numbering  at  a  node 
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node  ir  \,f  jr  common  math  language 


0000  V  at  lower  left  corner 
0  0  0  1  V.f 

0010  Vr 

0  0  !  1  Vrtf 

0100  V  at  lower  right  corner 

0101  V* 

0110  Vr 

0  I  1  1  Vry 

10  0  0  V  at  uppe'*  left  corner 

10  0  1  Vc, 

10  10  Vr 

10  11  Vr? 

1  1  00  V  at  upper  right  corner 

1101  V, 

1110  Vr 

1111  Vrf 

Table  3.4,  Numbering  of  nodes  and 
parameters  in  a  simple  quad. 

This  all  is  rather  irrelevant  to  the  casual  reader.  It  is,  however, 
essential  to  the  person  implementing  the  method  efficiently  on  a 
computer . 

Remark  •  The  coefficients  gu>  are  functions  of'-  rD  ,  rt  ,  y,  -  %.  It  follows 
that  the  coefficients  are  the  same  for  two  quads  sharing  these 
parameters.  If  all  simple  quads  in  a  'layer'  ra  (  r(  r,  have  the  same 
angular  width  -  fa,  then  the  coefficients  need  to  be  evaluated  only 


2  I 

3  I 

4  1 

5  2 

6  2 

7  2 

8  2 

9  3 

10  3 

1 1  3 

12  3 

13  4 

14  4 

15  4 

16  4 


once  1 
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3.4.2.  Partial  normals  of  a  compound  element 

Again  it  will  be  sufficient  to  consider  the  two~d i mens i ona I  case  and  to 
describe  the  procedure  on  hand  of  a  specific  example  It  is  the  compound 
element  described  above  in  section  3.1.6.  and  shown  in  figure  3.8 


Figure  3.8.  Compound  quad. 

We  start  by  interpolating  the  parameters  at  the  auxiliary  node  8  from 
those  at  nodes  4  and  5.  The  function  V(r,y)  is  given  in  the  upper  quad 
4,  5 ,  6,  7  by 

V(r,y)‘-  Z  £  (r-rt,  ?-*,)  C3.65) 

i.«*f  j-1 

It  is  useful  I  to  note  that  the  basis  functions  belonging  to  the  upper 
nodes,  i  .  e .  S"6j ,  <r7j  vanish  together  with  their  first  derivatives  at  the 
line  from  4  to  5.  Hence,  as  long  as  we  are  only  interested  in  zero-th 
and  first  derivatives  of  V(r,y)  at  this  line,  we  may  use  the  simplified 


.  iiit&.'ti, 


-  SQ 


reoresentat 


Vir.y)*  Z  Z  VL;  <u;  (r-rf,  -fj  +  0((r-rj2)  (3  66) 

j-i  J  J 

Using  a  more  explicit  represented  i on  we  may  write 


V(r.<j):  Z  [  Kw*1fy,(r'r,;r*'c^T.i,(»'& V*)  ♦ 

•  -  0 


(3.67) 


"  T1;  r  ( r-r, ;  >i  -  n )  •  (<f  -  y,;  ft-y,/ 


The  functions  o)r,;(x;h)  are  defined  in  section  3.  1. 2.,  equation  (3.8).  We 
are  now  able  to  interpolate  the  nodal  parameters  of  the  aux i ! ,crv 
node  8 


"  ^V.y  yjrVr  (  0  i  r'-r')rjH;h(  W  fz-%)  . 


(3.68) 


for  the  functions  >"Lltx)  see  equations  (3  9),  (3.18).  Switching  to  the 
more  condensed  no  tat  on 


i  -  2 1  +  i  +  ] 

i '  -n  2  '  +  '  x  • 

j  c-.<-  -y 


we  write  this  relationship 


(3.69) 


1/  . ,  3  JL  '  +  V- ;r;  V^. ;  ) 

?  )  I  .  J  J  ’J  J  J 

jC  1 


(3  70 
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In  an  even  more  condensed  form  we  write 

I 


v,  -  (W,  <4) 


K 

W 


V 


(3.71) 


where  U  is  an  8x4  matrix. 

The  next  step  is  to  form  according  to  section  3.4.1.  dealing  with 
simple  quads  the  partial  normals  for  the  lower  left  subelement.  Let 
these  equations  be  written  as 


*  tL 

^  i-L 

„  LL 

LL 

r 

|G„ 

&v. 

Q 

j 

LL 

Gz" 

r  u 
Gzt 

• 

G  UL 
<** 

G  tt 

0," 

I 

GS 

&LL 

PL 

CL 

1 

i 

[ 

K  ! 


=  o 


One  eliminates  the  auxiliary  node  by  the  substitution 


V  K  us) 


Ik 


VL 


(3  72) 


(3  73) 


Similar  to  a  parameter  transformation  in  adjustment  by  variation  of 
parameters,  we  must  consider  the  full  transformation  matrix  implied  by 

Vt] 

K 

t  i  i  v„ 


Vi 

Vu 


N  L  u,  * 


(3.74) 


-  6' 


V 


A  v 


and  transform  the  cdovo  normals  C3  72),  now  written  as 


G  V  --  0 


(3.75) 


according  to 


ArG  A  V  «  0 


(3  75) 


”he  resu It  :  s 


I 
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—  LL  LL  LL 

&V,  1  GVi  *  Gjy  U, 

Gi«  s  G#  ^4 


G  a  * 

r;  “•  y. 

GJV  *  U, 

r  u 

Gw(f  U 

:  t  4 

rr  ai 

r  a 
&2S  ' 

In  a  s i m i 

1  ar  way 

we  obta 

in  for 

the  1  ( 

! 

r  -  lr 

GIX 

—  LR 

G*, 

r  tR 

—  LR 
°IS 

-  LR 

—  LR 

G» 

—  LR 

Gu 

l 

-  LH 

0,x 

r  * 

g,;a 

—  LR 
GHS 

r 

?  LR 

G« 

r 

r  M 
Gjs 

By  a  d i rect  app 1 i 
upper  element 

cation 

of  the 

s  i  mp  1 1 

c 

c 

u 

o* 

c 

GSS 

O' 

Gr6 

G^ 

G* 

oj 

L/ 

Gis 

gJ 

g7^ 

G7/ 

r„  u.  ,  ,  r_  ll. 


(3.73) 


(3.79) 


(3  80) 


Now  the  partial  normals  (3.77),  (3.79),  (3.80)  of  the  3  subelements  are 
added  into  the  normals  of  the  compound  quads 
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r. 

\  ''n 

On 

Gic, 

G17  1 

v, ; 

j 

G1Z 

G2ti 

Gis 

G11  j 

!4; 

1 

G* 

r, 

Gy. 

G* 

G* 

0J> 

:  v,  i 

1  Gin 

0„ 

G^ 

G,s 

G« 

G\,>  j 

I 

•  1  K  i  C3 

81) 

o„ 

0n 

G„ 

G$s 

Gk 

Gn 

;  Yn 

* 

! 

C\z 

6*1 

G„ 

G«. 

!  vi 

i  1 

t - 

0 

"J 

G^ 

Gv 

(j-H 

G->s 

G» 

~  1 
G?7 

J 

VJ 

L  J 

Any 

of  the 

submatr 

i  ces  3 

Lj  ■  S 

c  sum 

V 

—  LL 

Gu  - 

&..tR 

LJ 

n  6 

V 

k 

C3 

82) 

where  Gk 

=  3,  GG*  =  3, 

O  is 

6  is  assumed. 

1  f  one 

of  the  nodes  i 

or  j 

does  not  par t i c 1  pa 

t  e  i  r 

the  subauad  cL 

_D 

or  'J . 

3.5 

More 

aener a i 

bcs  i  s 

fund. 

a  ns . 

This  section  Tay  se^ve  as  as  outl'ne  for  the  const~uction  of  shape 
functions  which  (!)  have  a  mere  general  analytical  r eoresentat i on  than 
piecewise  n-cubic  po  I  vnorn  •  a  !  s,  (2)  are  C  continuous,  and  (3)  are 
computat i anal 'y  si i ! •  reasonably  nff : c  e-t . 

The  not  i  vat ;  on  to  look  for  o  i  t  emat  ve  shape  functions  may  come 
from  the  des-e  to  use  funct'ons  ~h:ch  conform  to  the  potentials 
attenuation  ,n  a  better  wov .  it  appears  tempt  1 ng  to  replace  polynomials 

in  r  by  polynomials  m  !/r,  for  exam,: '  a  Th  1  s  ‘ooks  oar  1 1  cu  I  ar  I  y 

attract  vo  tor  elements  fc’/'er  away  from  the  earth  surface  It  also 

ooer>s  a  way  to  represent  '.he  or,  t  d t  a  i  -n  the  entire  ex  ter  •  or  of  some 


sphere  by  a  f  >  n  1  te  numbe'- 


' t . n 1 t e  size  Thi 


s  aspect  mi 


j  1  o"on 
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be  treated-  further  in  the  next  section. 

The  computat i ona I  efficiency  of  cubic  polynomials  comes  from  2 
sources.  Cl)  Such  polynomials  are  easily  evaluated,  (2)  the  integrals  to 
be  evaluated  in  order  to  form  the  partial  normals  of  a  simple  element 
decompose  into  a  moderate  number  of  products  of  integrals  over  a  single 
variable.  The  second  feature  extends  to  the  basis  functions  to  be 
outlined  below,  the  first  one  extends  only  partly. 

As  i  n  sect  ion  3.1.1.,,  one  starts  ago  in  with  the  unit  interval.  One 
introduces  functions  -yuCx)  such  that  the  j-th  derivative  of  Y'tj  Cx)  at 
the  i - tn- i n ter va I  end  is  unity,  while  all  other  derivatives  at  this  and 
the  other  interval  end  vanish.  It  is  important  to  note  that  the  highest 
derivative  considered  this  way  need  not  be  the  same  at  both  interval 
ends.  To  have  something  specific  in  mind,  consider  the  basis  functions 

y  (y)  *  3 x1'  -  h  x-3  +  1 

J  00 

yo1  (*)  -  +x 

yoz(x)  =  *  z*2'  (3. 83) 

Yio  W 

r„  M 

In  this  example  the  second  derivative  is  only  modelled  at  the  left  end 
of  the  interval,  but  not  at  the  right  end.  If  the  above  basis  functions 
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are  changed  to 


'Y'  (x)  -■  -  6  xs  +  15 -  10 

l  OO 

Yo  (x)  *  -  3xr  +  Sx1*  -  6  x1  + 

ruW--K*Tx'- tv3+tx" 

Yo  (x)  *  6xs- l^x1* +10x2 
Y"  (*)  -  -  3xs  +  7  -  kx1 


(3.8^) 


■'  T 


d  S’ 

‘  v  - 


V  -1  3 

X  +  —X 

z 


then  we  are  dealing  with  Herrnile  quint ic  polynomials.  At  each  node  the 
zero,  first  and  second  derivative  can  be  prescribed. 

With  the  chosen  set  of  basis  functions  for  the  unit  interval  one 
can  proceed  to  any  finite  interval  by  (cf.  section  3.1.2): 


T  (V;h)  ■  h'y. .(£) 


(3.85) 


As  in  section  3. I. 3.,  two  dimensional  basis  functions  are  again  obtained 
by  forming  products 


T L,iyiJy{Ky) '  rij(y;a^Tij(y;b)  C3$6> 

As  indicated  by  the  bar  over  the  last  y,  there  is  no  need  to  have  the 
same  type  of  basis  functions  for  x  and  for  y.  For  example,  one  may  use 
quint ics  in  x  and  cubics  in  y. 

In  order  to  ensure  C"1  continuity,  the  type  of  basis  functions  used 
in  the  x-direction  should  be  matched  by  the  element  adjacent  in  the  y 
direction  and  vice  versa.  This  poses  some  problems  if  one  fuses  simple 
elements  to  compound  elements  However,  as  long  as  the  basis  functions 
are  polynomials,  there  is  normal ly  a  way  out  of  any  di lemma.  In  any 
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case,  (^continuity  must  be  checked  very  carefully  in  ail  unusual 
s i tuat i ons . 

Remark ■  If  quintics  are  used  consistently,  the  resulting  function  is 
even  C1"  continuous.  This  may  look  attractive.  Among  others,  it  offers  a 
way  to  force  the  Laplacean  to  be  zero  at  al i  nodes.  Houever,  a  heavy 
price  has  to  be  paid  in  terms  of  d1  parameters  per  node  (d  .  .  .  number  of 
dimensions).  For  d  =  3  there  are  27  parameters  per  node,  while  in  case 
of  bi cub ics  there  are  only  8. 


The  nodal  parameters  needed  in  these  more  general  finite  element 
representations  follow  automatically  from  the  one-dimensional  basis 
functions  used.  Two  dimensions  are  typical  enough  for  our  outline  here 
Any  node  has  a  degree  in  the  x  and  y-di recti  on.  These  two  degrees  need 
not  be  equal,  but  they  pose  restrictions  to  the  basis  functions  of 
adjacent  elements.  Suppose  the  degree  is  2  in  the  x-direction  and  I  in 
the  y~di recti  on.  This  holds  if  quintics  in  x  are  multiplied  with  cub ics 
in  y.  Then  the  2x3  nodal  parameters  are  automatically  given  by  all 
symbolic  products  of  any  two  elements  out  of  the  sets 


(3.87) 


i  . e .  they  are : 


(3.88) 


The  integration  procedure  outlined  in  section  3.4. f.  and  the  partial 
normals  of  a  simple  element  must  be  modified  if  other  basis  functions 
than  cubic  polynomials  are  used.  The  modifications  are,  however, 
moderate.  They  amount  to  a  change  of  the  integration  formulas  for  the 
required  products  of  one-dimensional  basis  functions  and  their 


Figure  3.9.  Sample  basis  functions  for  an  interval 
extending  to  infinity 


We  combine  them  e.g.  with  cubics  in  y> ,  which  we  denote  by 


Tij  &%>%'%)  %  6  y  *  %  (3  91) 

Hence  we  deal  with  the  basis  functions 

tjf  *  0  j  j*  ,  Ly.  jy  ~  0,  1 

There  is  no  C  continuity  problem  if  an  element  partition  is  chosen  as 
shown  in  figure  3.10  and  if  bicubics  are  used  for  all  elements  of  finite 
size. 


Infinite  quad 


Figure  3.10.  Element  partition  making  use  of 
infinite  el ements 
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Next  we  outline  the  integration  procedure  in  detail  for  the  infinite 
element.  Confer  the  parallel  developments  of  section  3.4.1  on  the 
partial  normals  of  a  simple  quad .  The  representat i on  of  the  potential 
V(r,yO  in  air,  y  <  is 


V(ry)-  Z  VLryJrV  X£r.Jr.  (r)  Thy  {<?%;%■%)  (3 


S’jVjV 


93) 


The  index  i '  is  fixed  at  zero  (we  retain  it  only  for  ease  of  comparison 
with  section  3.4.1.,  equation  (3.38)).  All  other  indices  run  from  0  to 
1 .  The  Laplacean  is: 

Av(r^  ’ .  K,  ,Wv  { ’fiWr- (r)  r.vjy  <»-»■•**»■•> + 


<-r  ‘■'/Jr  J<f 

*TX,irVlr)r.w 
*  m  Tiijiy  (?-*••■ 

This  equation  is  abbreviated  as 

5  M  hw> 
’*yW  JT, 1  VJr  W 


(3.94) 


AV(ry)  •  £  v,.  Zf'*',.(r)  h‘k''  (y) 

rUfJrJ,f  t,‘Jr  i»jf 


(3  95) 


The  principle 
«>  y. 


1 :f  f  (A'/(r</))Vd.rdy  --  Min 

r-r0 

leads  to  the  partial  normals  of  the  infinite  quad: 


(3.96) 


To.  .  ■  .  .  .  V  •  .  . 

JrJy /  •-r'l-y'jr'iy1 


0 


(3.97) 
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In  analogy  to  section  3.4.1  ,  the  coefficients  follow  from  the  formula- 

3 

_p  KK  -f-KK 

(3.98) 


■‘"flrJ'j 

The  integral s 
kk' 


-j-  KK  -j-  KK 

k.k'.|  LrJr  lrV  WWJ* 


r . 

<.y)y  uy  Jy* 

are  unchanged.  On  the  other  hand: 

Cw  '  /<^<«iv/r(r,X^MlrV(r'  r 


e(W*e(K*M 


(3.99) 


Table  3.3  of  section  3.4.1.  for  the  dr,  e  is  still  valid.  If  one  writes 

the  Xj.ri.rJr  '  n  the  form 

1  r  2.  +  jr+  2 

i*Vrv(r)-  5  T^tt  <3;005 


then  a  compact  3xlx2x2  table  of  the  coefficients 

o  i  aLr  i  Z 

r  .  *  Cr  1  0 

<J-rLrJrC  o  4  j>  if 

Oit  *  4 


(3. 101) 


may  be  stored  in  the  computer.  This  CHC-7ob  s  the  counterpart  to  the 
PSC-Tab!e  introduced  in  section  3.1.1.  help  of  the  CHC-Table 

one  computes: 


■^i-rjr^yiy  ^  Ccir(K)  lr  Jr  «  C  dr(K'J  Lr' Jr' ^ 


(3.102) 


jr  fjr'-  c£r(K)  -  o(r(K')r 


pKt*'Z'~  1 
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Thereby 

PKgK.6‘-  =■  1*  oir CK)  +  i  *  1+  d.r(K‘)  +  e'~  e(K)~  e(K‘)-1  (3.103) 


3.7.  The  outer  zone. 

The  infinite  outer  space  of  the  earth  can  not  be  partitioned  into 
infinitely  many  elements  of  finite  size.  The  field  is  most  detailed  in 
the  vicinity  of  the  surface.  Hence  fairly  small  elements  must  be  chosen 
there.  The  elements  may  increase  as  we  go  outward.  The  increase  in  size 
is  partly  achieved  by  the  natural  increase  of  the  surface  element  rdrdy 
or  r1cosydr  dycU  in  polar  coordinates.  To  a  greater  extent  it  is 
achieved  by  a  lumping  of  elements  as  discussed  in  the  section  on 
compound  elements.  At  a  certain  level  the  finite  element  partition  may 
stop  at  all.  The  field  is  then  represented  in  a  spherical  shell 
(circular  ring).  It  is  forced  to  be  (nearly)  harmonic  there  by 
minimizing  the  integral  over  squared  Laplacean.  It  is  now  necessary  to 
ensure  in  some  way  that  the  field  in  the  outer  space  of  the  shell  is 

(1)  consistent  with  the  field  in  the  interior  of  the  shell, 

(2)  approaching  zero  if  the  radius  tends  to  zero,  and  finally 

(3)  harmonic. 

Consistency  means  in  a  strict  sense,  that  the  outer  field  is  the 
analytical  continuation  of  the  field  within  the  shell.  In  a  finite 
element  context  this  requirement  must  be  somewhat  relaxed.  Not  even 
within  the  shell  do  we  have  a  strictly  analytical  field.  However,  if  we 
have  C"*  continuity  there,  as  we  do  when  cubic  polynomials  are  used, 
approximate  C1  continuity  across  the  outer  boundary  of  the  shell  is  a 
minimum  requirement. 

There  are  various  ways  to  deal  with  the  problem  of  the  outer 
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field.  We  sha I  I  out  I  i ne  them  now , 

Cl)  Choose  the  radius  of  the  outer  shell  large  enough  and  force  the 
field  and  its  first  derivatives  to  be  zero  there.  A  large  enough  radius 
means  that  the  size  of  field  and  its  derivatives  must  be  below  the 
desired  computational  accuracy.  It  must  be  borne  in  mind  that  our  field 
is  actually  a  disturbing  field  superimposed  upon  a  reference  field. 

(2)  Partition  of  the  outer  space  into  finitely  many  elements  of  infinite 
size.  Refer  to  the  previous  section  3.6  on  elements  extending  to 
infinity  in  one  direction. 

(3)  Use  of  Greens  formula.  If  potential  and  (outer)  normal  derivative 
are  prescribed  at  the  boundary  P,  then  the  potential  in  the  exterior  of 
r  can  be  represented  by  Greens  3r<i  identity.  CSee  e.g.  He  i  skanen-Mor  i  tz 
(1967),  equation  1-29'). 

^yr^(P)  -  j'  {  V(Q)  ~s^(  <(Piq))  dT,  Jl r>,a)  J  dria)  (3,1045 

r 

This  formula  refers  to  the  3-dimensional  case,  I  CP , Q )  is  the  distance 
between  the  reference  point  P  and  the  point  of  integration  G.  Confer 
f i gure  3.11 


Figure  3.11.  Explains  the  definition  of 
outer  normal  n  and  distance  I 
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Consider  now  an  element  partition  within  a  shell  around  the  earth 
surface  as  shown  in  figure  3.12. 


Figure  3.12.  Fusing  the  inner  and  the  outer  zone 

The  surface  P  of  integration  in  equ.  (3.104)  is  assumed  spherical 
(circular)  and  situated  slightly  inside  the  outer  boundary  of  the 
shell.  Greens  third  identity  may  now  be  applied  to  any  node  P  situated 
at  the  outer  boundary  denoted  P„.  Since  V C Q )  and  is  expressed 

linearly  in  terms  of  the  parameters  belonging  to  nodes  situated  in  the 
layer  between  P0  and  Pt ,  we  get  linear  equations  relating  V(P)  to  nodal 
parameters  of  the  outer  layer  (bounded  by  PL  and  Pa) .  Let  V(Q)  be 
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represented  as 


v(a)  •  Z  Z  vu  SJQ) 

l*l  j  y  y 


(3.105) 


Uhere  L  represents  the  set  of  nodes  par t i c i pat i ng  in  the  outer  layer, 
and  j  runs  over  the  parameters  of  any  node.  The  Sij(Q)  are  the  shape 
functions  introduced  in  section  3.2.  We  get 


3V(Q) 

dn 


Y  V  V  jMgj 

ui  j  ^ 


(3. 106) 


vw-iz{{  Z  ]dP(s> 


The  amounced  linear  equations  are  thus  obtained; 

v(p)  *  Z  Z  c  V 

16L  J  J  J 


(3.108) 


In  a  similar  manner  we  may  even  express  derivatives  of  V  at  P. 

D i f ferent i at i on  of  Greens  formula  poses  no  problem  because  P  is  not 
situated  at  P.  Therefore  l(P,0)  never  becomes  zero  For  example  we  may 


eva  i  u'-  ;e 


dv(P) 


1  Z  a.;y: 


5viP)  sr  \  / 

<?A  *  lr 


(3.109) 


dxv{p) 


Representing  V(P)  and  the  ’horizontal  derivatives’  of  V  ot  P  in  this 
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way,  and  doing  this  for  all  nodes  of  PG,  provides  us  with  a  set  of 
linear  equations  which,  together  with  the  normal  equations  of  the 
interior  of  the  shell,  will  ensure  a  potential  which  is  "reasonably 
harmonic"  outside  the  earth.  A  disadvantage  of  this  method,  which  is 
proposed  in  McDona I d-Wex I er  (1972)  and  further  discussed  in  Zienkiewicz 
et  al .  (1979),  is  the  need  for  the  evaluation  of  a  number  of  integrals 
which  are  not  ali  too  simple.  The  following  3  remarks  are  considered 
important . 

Remark  I ■  The  reader  may  wonder  why  we  used  only  V  together  with  its 
horizontal  derivatives  in  the  above  compatibility  equations.  The  answer 
consists  of  two  parts.  First  we  have  to  point  out  that  the  underlying 
assumption  is  that  of  a  tricubic  representat i on  of  the  field.  The 
"horizontal"  nodal  parameters  V,  Vy,  V* ,  are,  so  to  speak, 

responsible  for  the  variation  of  V  along  the  surface  P0.  Secondly,  we 
should  stress  that  our  finite  set  of  equations  substitutes  for  an 
infinite  set  of  equations  (3.107)  in  which  P  varies  all  over  P0.  The 
variation  of  P  and  V(P)  along  P0  is  now  logically  replaced  by  the 
variation  of  those  parameters  in  the  cubic  r epresentat i on  which  are 
responsible  for  the  behaviour  of  VCP),  for  PePa. 

Remark  2-  The  discussion  in  Remark  I  also  demonstrates  an  alternative 
way  in  which  Greens  third  identity  can  be  used  in  order  to  ensure 
compatibility  of  the  field  across  P„.  Instead  of  using  equations  (3.109) 
for  Vy,  V*,  Vy>  we  could  just  use  equation  (3.108)  for  V(P);  however,  we 
must  o  it  4  times  as  often  as  there  are  nodes  on  P0.  A  way  to  do  this 

would  to  use  (3.108)  for  the  nodes  P  on  P„  and  also  for  the  points 
halfway  between  two  adjacent  nodes.  Confer  figure  3.13. 
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x  x  x  x  x 
0X0X0 
X  X  X  X  X 
0X0X0 

finite  element  nodes  on  P8 
additional  points  p  at  which 
Greens  identity  is  evaluated 

Figure  3.13.  Pattern  of  nodes  as 
explained  in  the  text 

Remark  3:  One  could  even  go  one  step  further.  One  could  use  Greens 
identity  for  more  horizontal  parameters  than  this  is  implied  by  the 
discussion  in  the  foregoing  remarks.  Exact  fulfillment  of  the  identities 
then  can  no  longer  be  postulated.  One  would  apply  weights  to  the 
discrepancies  and  add  their  weighted  sum  of  squares  to  the  other 
functional  to  be  minimized  In  this  way  one  even  ends  up  with  a  positive 
definite  symmetric  system. 


o 

x 


(4)  Combination  with  spherical  harmonics.  Imagine  the  potential 
represented  at  the  outer  boundary  P0  in  terms  of  spherical  harmonics 


v(p)-  t-rkr'i  Cnm  Hnm(P) 

niO  r  m»-n 


(3110) 


We  use  the  symbol  "H*  because  the  “S*  is  already  reserved  for  the  shape 
functions.  We  may  also  represent  the  horizontal  derivatives  of  V  in  this 
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way 
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(3.111) 


We  form  these  equations  for  any  node  at  ("V  Imagine  that  there  are  just 
as  many  equations  as  there  are  coefficients  Cn»,.  Then,  hopefully,  the 
Com  could  be  evaluated  from  them  However,  we  do  not  propose  to  do  this 
Instead  we  propose  to  add  equations  for  the  remaining  nodal  parameters, 
i  .  e . 
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In  this  way,  we  will  ensure  that  the  potential  is  reasonably  harmonic 
outside  of  P„  compatible  with  the  field  within  the  shell,  and  tending  to 
zero  as  t — . 

Remarks  I  to  3  given  above  under  (3)  apply  mutatis  mutandis  to  the 
present  situation. 


3.8.  Local  data  deficiencies 


It  may  happen  that  data  are  abundant  and  redundant  in  some  areas  while 
in  others  they  are  sparse  and  deficient.  Lack  of  sufficient  data  in 
local  areas  can  cause  rank  deficiencies  or  near  singularities  of  the 
normal  equations.  To  illustrate  our  point,  consider  Dirichlet's  problem 
for  the  exterior  of  the  unit  circle-' 

AV  (r.  <f)  -  0  ;  r  >  1 

V  (1.<f)  ■  f((f)  (3. 113) 

V(r,<f)  •  0( Log  r)  ;  r  -*  <*> 

Suppose  that  element  partition  and  distribution  of  data  are  as  shown  in 
f i gure  3.14. 


Figure  3.14.  Distribution  of  data  leading  to 
singularities  of  the  normal  equations. 

Infinite  elements  are  used  for  the  outer  zone.  It  is  seen  that  data  are 
mi.ssing  or  def >  int  in  3  of  the  8  intervals  at  the  c i rcur»ference  of  the 
unit  circle.  It  follows  that  a  cubic  spline  function  s(y)  approximating 


i 


J 
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fCy)  in  the  least  squares  sense  is  undetermined.  Hence  the  contribution 
of  the  measurements  to  the  normal  equations  will  have  a  rank  smaller 
than  in  the  case  of  sufficient  data.  The  addition  of  the  contribution  of 
the  field  to  the  normal  equations  does  not  remedy  the  situation.  The 
field  contribution  merely  serves  to  ensure  (approximate)  harmomcity  in 
the  exterior  region  r  >  1 .  If  the  boundary  data  are  undetermined  in 
Dirichlet's  problem,  the  whole  field  will  be  undetermined. 


Remark •  Only  to  freshmen  of  adjustment  courses  it  may  appear  paradoxical 
that  a  rank  deficiency  occurs  when  the  number  of  observations  exceeds 
the  number  of  unknowns.  Think  of  a  network  which  may  be  rigid  and  stable 
in  some  areas  and  poorly  determined  in  others 


There  are  at  least  two  ways  to  remedy  the  rank  deficiencies  due  to 
lack  of  data  in  some  areas,  (t)  The  elements  may  be  chosen  larger  there. 
(2)  A  third  contribution  to  the  normals  may  be  calculated  which  results 
from  the  square  of  a  certain  curvature  norm  applied  to  the  spline  s(y) 
approximating  f(<p). 

Approach  (1)  is  not  too  strongly  recommended.  It  may  cause  a  loss 
of  regularity  in  the  element  partition.  As  a  consequence  the  computer 
programs  could  become  more  involved. 

Approach  (2)  is  related  to  least  squares  collocation  using 
splines.  We  take  the  finite  element  representation  of  V(1,y) 

Vtt.y)--  ZI  Z  Vii  Stj il.y)  (3.M4) 

The  index  set  L  comprises  the  nodes  situated  on  the  circumference  of  the 
unit  circle.  We  take  the  second  derivative 


d*V(i.<r) 

dy1 


’  z 


SCi.y) 


(3. 1 15) 
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and  consider  the  integral 
•Zjt 

-J-  E  ’  f  [Z  E  Vy  ^  Stt.y)  %(v>oly  0.116) 

Here  p(y)  is  a  weight  function  which  can  be  taken  p(^f)  =  I .  It  may  also 
be  taken  positive  in  areas  of  poor  data  and  zero  in  areas  of  sufficient 
data.  Taking  the  variation  of  (3.116)  with  respect  to  Vy  yields,  very 
much  in  the  same  way  as  outlined  in  section  3.3.2.,  a  third  contribution 
to  the  normals.  Sparsity  of  the  normals  will  not  be  impaired  in  any  way. 

The  method  readily  generalizes  to  3  dimensions.  The  second 
derivative  with  respect  to  y  may  be  replaced  by  the  surface  Laplacean 
thereby . 


Our  normal  equations  can  be  grouped  according  to  nodes.  If  d  denotes  the 
dimensions,  i.e.  d  =  I,  2  or  3  and  if  d-cubics  are  used,  then  for  each 
node  there  are  2d  rows  and  2d  columns  corresponding  to  the  2^  para¬ 
meters  per  node.  During  any  step  in  reordering  or  reducing  the  normal 
equations  the  2d  equations  for  one  node  will  always  be  lumped  together. 
We  may  view  the  system  as  composed  of  2dx2d  submatrices  which  could  be 
called  generalized  coefficients.  Elimination  is  then  carried  out  using 
these  generalized  coefficients  instead  of  conventional  scalars. 

In  the  original  normals,  two  nodes  are  coupled  by  nonzero 
off-diagonal  coefficients  if  there  is  a  finite  element  on  whose  boundary 
both  nodes  are  located.  Hence  the  system  is  sparse.  During  elimination 
the  coupling  increases  due  to  fill-in.  If  a  node  is  eliminated  all  its 
neighboring  nodes  which  are  not  vet  el iminated  become  coupled.  The  whole 
idea  of  sparse  elimination  is  to  renumber  the  nodes  in  a  way  that 
fill-in  is  effectively  kept  down.  We  shall  use  a  technique  which  is  a 
combination  of  nested  dissection  and  Helmert  blocking.  The  philosophy 
behind  this  procedure  was  discussed  in  great  detail  in  Meissl  (1980), 
sections  3.5.4  -5.  Here  we  shall  outline  it  on  hand  of  a  two  dimensional 
problem  with  element  partition  shown  in  figure  4.1. 


Figure  4.1.  Element  partition  serving  the  outline  of 
nested  dissection  in  conjunction  with  Helmert  blocking 
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Mind  that  we  work  in  the  plane  polar  coordinate  system.  Hence  the  region 
of  figure  41  is  actually  a  circular  ring,  the  left  and  right  boundary 
line,  and  the  nodes  located  on  them  are  to  be  identified.  We  first 
describe  the  block  design  in  a  bottom  up  fashion.  At  the  first  stage,  we 
consider  8  blocks.  Each  one  is  formed  out  of  2  adjacent  squares  located 
at  the  bottom  level.  Such  a  block  looks  as  shown  in  an  enlarged  way  in 
figure  4.2.  We  imagine  the  normals  formed  for  the  two  simple  quads 
composing  the  block,  and  the  normals  added.  We  eliminate  from  them  the 
node  indicted  by  a  circle  in  figure  4.2. 


Fig.  4.2.  First  stage  block 

We  call  this  the  inner  node  of  this  block.  The  other  nodes  C indicated  by 
crosses)  are  called  junction  nodes.  We  obtain  a  partially  reduced  system 
of  5  junction  nodes. 

At  the  second  stage  we  enlarge  the  blocks  of  stage  one  by  adding 
two  adjacent  squares  of  the  next  higher  layer.  The  block  looks  as  shown 
in  figure  4.3.  There  are  still  8  blocks. 


Figure  4.3.  Second  stage  block 
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Each  block  again  shows  one  inner  node.  However,  this  time  there  are  7 
junction  nodes.  The  normal  equations  for  such  a  node  are  obtained  by 
adding  the  partially  reduced  normals  of  stage  I  and  the  normals  for  the 
upper  squares.  After  elimination  of  the  inner  node,  we  obtain  a 
partially  reduced  system  for  the  7  junction  nodes. 

At  stage  3  we  adjoin  to  the  block  of  the  previous  stage  one 
compound  element  of  the  third  layer.  Ue  obtain  8  blocks  as  shown  in 
f i gure  4.4. 


* - ¥ 

:c  :: 

— - O - 

:c  :: 

x  * 

Figure  4.4.  Third  stage  block 

The  normals  are  formed  by  adding  the  partially  reduced  normals  of  stage 
2  and  the  normals  for  the  compound  upper  quad.  The  inner  node  is 
el iminated. 

At  stage  4  we  join  two  adjacent  blocks  of  the  previous  stage  and 
eliminate  the  4  inner  nodes.  Confer  figure  4.5.  The  number  of  blocks  is 
4  in  this  stage. 
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Figure  4.5.  Fourth  stage  block 


At  stage  5  we  adjoin  to  the  previous  block  the  compound  elements  of  the 
upper  layer.  We  obtain  blocks  as  shown  in  figure  4.6.  The  scale  of  this 
figure  is  now  the  same  as  that  of  figure  4.1 .  The  number  of  blocks  is 
sti I  I  4. 


Figure  4.6.  Fifth  stage  block 


At  stage  6  we  lump  two  adjacent  blocks  of  stage  5.  See  figure  4.7. 
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Figure  4.7.  Sixth  stage  block 
At  stage  7  we  do  a  similar  thing;  confer  figure  4.8 


* 


Figure  4.8.  Seventh  stage  block 


Recall  that  the  nodes  of  the  two  outer  boundaries  are  actually 
identical.  In  the  final  and  8-th  stage  we  eliminate  just  these  nodes 


Figure  4.9.  Eights  stage  block 


The  Helmert  blocking  procedure  is,  as  we  feel,  best  understood,  if  it  is 
described  in  the  above  bottom  up  fashion.  In  order  to  analyse  various 
block  designs,  however,  in  particular  block  design  for  3-dimensional 
problems,  we  prefer  a  top  down  fashion,  which  we  briefly  outline  for  the 
above  example. 

At  stage  8  we  slice  the  spherical  ring  by  a  dissecting  fine.  Inner 
nodes  are  those  located  at  this  line.  The  dissecting  line  is  indicated 
by  an  arrow  in  figure  4.9. 

At  stage  7  we  dissect  the  block  again  (arrow  in  figure  4.8)  and  we 
repeat  this  at  stage  6.  Inner  nodes  are  still  those  nodes  located  at  the 
dissecting  line.  At  stage  5  the  dissecting  line  is  horizontal  (figure 
4.6).  Here  we  discover  the  rule  which  will  be  general  enough  for  all  our 
e I ement  par  t i t i ons ■ 

Rule  for  identifying  inner  nodes  and  junction  nodes :  Inner  nodes 
are  those  nodes  located  at  the  dissecting  line  which  have  not  been  inner 
nodes  before.  (In  3  dimensions  we  deal  uith  dissecting  surfaces  rather 
than  lines).  Junction  nodes  ore  inner  nodes  of  the  next  higher  stage 
together  with  junction  nodes  of  the  next  higher  stage  located  in  the 
block  under  consideration. 

The  dissecting  lines  of  the  stages  1-4  are  easily  recognized  from 
figures  4.2  through  4.5  (arrows). 
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4.2.  Global  soiut i on 


We  consider  a  finite  element  partition  which  allows  to  model  the  details 
of  the  field  near  the  surface  with  about  the  same  resolution  as  a 
conventional  surface  layer  or  collocation  solution  based  on  1° x  1°  block 
averages  of  gravity  anomalies. 

The  dimension  is  3  and  we  use  tricubics.  This  implies  that  there 
are  4  parameters  available  at  each  node  which  are  responsible  for  the 
hor i zon ta I  variation  of  V  in  the  vicinity  of  this  node.  Hence  the 
element  size  near  the  surface  must  be  chosen  as  2° x  2° .  Mind  that  the 
number  of  blocks  is  the  same  as  the  number  of  nodes.  Consequently  there 
are  4  parameters  available  to  model  the  horizontal  variation  of  V  in  a 
block.  It  follows  that  the  averages  over  gravity  anomalies  in  the  four 
1° x  1°  subblocks  can  be  matched  exactly.  Our  surface  of  computation  will 
be  a  sphere.  It  may  be  imagined  as  a  sphere  slightly  below  the  earths 
surface.  The  computational  effort  of  solving  the  normals  is  unaffected 
by  the  use  of  a  more  complicated  reference  surface.  Only  the  formation 
of  the  normals  takes  longer.  Asymptotically,  i.e.  for  very  large 
systems,  formation  is  negligible  in  comparison  to  solution  By  assuming 
a  concentric  sphere  with  an  appropriately  larger  radius  we  obtain  a 
spherical  shell.  Our  element  partition  refers  to  this  shell.  It  is 
specified  in  figure  4  10(a), (b)  and  in  tables  4,1  to  4.3. 


zone 

from  to 
(degr .  I  at i t .  ) 


block  size 
I  at i t .  ( ong i t . 

(degrees) 


Number  of 
b  I  ocks  a  I ong 
a  para  I  I e I 


Number  of 
blocks  in 
zone 


0  72  2  2  180  6480 

72  84  2  6  60  360 

84  90  6  18  20  60 


Table  4  1  Element  part ;t ion  in  first  (bottom) 
layer  as  wel !  as  in  second  layer 
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F i gu re  4 . 10(a) 

Element  par  till  on  for  global  solution 
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Figure  4.10(b).  Element  partition  for  global  solution. 

Detailed  view  of  area  near  the  pole 


The  tables  apply  to  the  northern  hemisphere.  The  partition  for  the 
southern  hemisphere  follows  by  symmetry.  In  the  bottom  layer  there  are 
13800  blocks  over  the  whole  sphere.  Note  that  the  area  of  the  sphere 
divided  by  the  area  of  a  2°x  2°  equatorial  block  gives  about  10300. 

Let  us  pause  for  a  moment  and  reflect  on  the  computational  effort 
associated  with  a  surface  layer  solution  based  on  2°x2°  blocks.  (Recall 
that  only  l°xl°  blocks  would  be  about  equivalent  to  our  finite  element 
treatment).  We  are  dealing  with  a  system  of  13800  unknowns.  The  unknowns 
are  the  densities  for  the  blocks.  To  solve  a  linear  symmetric  system  of 
i  unknowns  by  a  direct  elimination  method  requires  about 
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(4.1) 


6 

elementary  steps.  (Confer  appendix  A,  equation  (A.6)).  One  elementary 
step  comprises  a  multiplication  followed  by  an  addition.  We  arrive  at 

IM22I  Z  Ell  (4  2) 

<5 

steps.  Assuming  a  computer  that  can  perform  500,000  steps  per  second 
(OSU  has  such  a  computer  at  the  present  time),  we  estimate  a  total 
effort  of 

243  hours  CPU  lime  (4.3) 

If  the  blocks  are  chosen  1° x  1°  the  number  of  blocks  multiplies  by  4. 
Hence  the  estimates  of  (4.3)  multiplies  by  4 3  =  64  giving  about 

15,000  hours  CPU  time  (4.4) 

We  return  to  the  finite  element  partition.  The  experience  gained  from 
the  computer  experiments  documented  in  chapter  6  and  dealing  with  the 
two  dimensional  Stokes  problem,  persuades  us  to  assume  two  layers  of 
elements  of  the  size  shown  in  table  4.1.  At  the  third  layer  blocks  are 
fused,  mostly  9  into  one.  Table  4.2  shows  the  block  size  for  the  third 
layer.  (Cf.  again  figure  4. 10(a), (b)) 


zone 

from 

(degr . 

to 

1  at  i  t . 

block  size 

lat.  long. 

)  (degrees) 

No.  of  b 

along  a 
paral lei 

locks  No.  of  blocks 

in  zone 

0 

72 

6 

6 

60 

720 

72 

84 

6 

18 

20 

40 

84 

90 

6 

90 

4 

4 

Table 

4.2. 

Element  part i t i on 

in  third  layer 

-  91  - 

At  the  fourth  layer  the  same  partition  is  used.  At  the  fifth  layer  we 
fuse  again;  table  4.3  shows  how. 

zone  block  size  No .  of  blocks  No.  of  blocks 

from  to  lat.  long.  along  a  in  zone 

Cdegr .  I  at  it.)  (degrees)  parallel 

0  72  18  18  20  80 

72  84  12  90  4  4 

84  90  6  360  I  I 

Table  4.3.  Element  partition  in  fourth  (upper)  layer 

Within  the  main  zone  -72°  to  +72°  the  basic  building  block  is  a  18°x  18° 
configuration.  We  call  it  a  6m-stage  block.  The  meaning  of  this  name 
will  become  clear  later  on.  The  face  of  such  a  6'h-stage  block  shows  a 
pattern  of  nodes  as  depicted  in  figure  4.11.  We  count  46  nodes  on  this 
face,  and  8  along  a  side  line. 

In  the  following,  we  assume  a  blocking  strategy  for  the  formation 
and  solution  of  the  normals  as  outlined  in  the  previous  section.  The 
normals  result  from  the  contribution  of  the  field  and  of  the 
observations.  Confer  section  3.3.  It  is  important  that  the  observations 
are  1  oca  I .  Any  observation  must  only  involve  points  situated  in  one  and 
the  same  element.  The  contribution  of  the  field  in  the  outer  zone  is 
assumed  to  be  taken  into  account  by  finite  elements  of  infinite  size. 
Confer  section  3,6.  The  partition  is  that  one  shown  for  the  last  inner 
layer  in  table  4.3. 
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The  blocking  will  be  hierarchical.  Block  boundaries  will  always  be 
composed  of  element  boundaries.  We  proceed  to  describe  the  design  and  to 
estimate  the  computational  effort  in  a  top  down  manner,  starting  at  the 
last  stage. 

At  the  last  stage,  which  is  stage  13,  there  are  only  inner  nodes. 
We  assume  that  they  are  comprised  of  all  nodes  situated  in  the 
equatorial  plane.  See  figure  4.12(13).  Figure  4.11.  shows  the  pattern  of 
nodes  which  repeats  itself  20  times  along  the  equator.  Hence  the  number 
of  nodes  i s 

i  =  20*46  -  20*8  =  760 

The  correspond i ng  s/stem  of  partially  reduced  normal  equations  has  8*i 
equations  because  we  have  8  parameters  per  node.  In  agreement  with 
equation  (4.1)  the  computat ional  effort  in  steps  is 

a  1.1SE10 

6 

translated  into  CPU  time  (assuming  500,000  steps/sec)  this  gives 

CPU  time  =  20.8  hours  [stage  13]  (4.5) 

At  the  next  lower  stage,  i.e.  stage  12,  we  have  2  blocks,  the  northern 
and  the  southern  hemisphere.  We  cut  off  the  polar  caps  at  latitudes  *84 
degrees  (cf.  figure  4.12(12)).  We  see  that  for  one  block 

i  =  248 

whereas 

j  =  760 

from  the  previous  stage.  According  to  equation  (AS)  of  appendix  A  the 
computational  effort  to  eliminate  n  interior  equations  in  the  presence 
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of  m  junction  equations  is 


n 2  n  m  (n  r  m) 

6  2. 


(4.6) 


In  our  case  n  =  248*8,  m  =  760*8.  Hence  the  effort  in  steps  is 


O'  (ifj) 
1 


(4.7) 


Translated  into  computer  time,  and  multiplied  by  the  number  of  blocks 
(=2),  this  gives 


CPU  time  =  55.5  hours  [stage  123 


(4.8) 


We  proceed  to  the  next  lower  stage,  i.e.  stage  I  I .  We  are  still  dealing 
with  2  blocks.  Each  block  is  divided  into  two  parts  by  the  Greenuich 
meridian.  Keeping  figure  4.12(11)  in  mind  we  count  in  one  block 

i  =  8*46  +  2*34  -  12*8  =  340 


interior  nodes  (along  the  dividing  meridian)  and 

j  =  760  +  248  =  1008 

junction  nodes  (at  y  3  0°,  i84°).  The  CPU  for  both  blocks  is  estimated 
at 


CPU  time  =  135.1  hours  [stage  II]  (4.9) 

At  the  next  lower  stage,  which  is  stage  10,  we  are  dealing  with  4 
identical  blocks.  Cf .  figure  4.12(10).  Each  one  comprises  a  quarter  of 
the  sphere.  We  subdivide  each  block  by  the  central  meridian.  For  one 
such  block  we  have 


i  =  4*46  +  34  -  6*8  =  170 
j  =  10*46  +  2*(4*46  +  34)  +  2*70  -  22*8  3  860 


-  96  - 


We  obtain  a  total  CPU  for  this  stage  os 

CPU  time  =  86.6  hours  [stage  10]  (4.10) 

At  the  S-th  stage  we  have  8  blocks  corresponding  to  the  8  octants.  We 
cut  off  the  region  near  the  pole  by  dissecting  along  the  surface 
‘f  =  -72°.  See  figure  4.12(9).  We  have 

i  =  5*46  -  6*8  =  182 
j  =  13*46  +  2*34  +  70  -  16*8  =  608 

we  obtain 


CPU  time  =  101.8  hours  [stage  9]  (4.10a) 

From  now  on  we  describe  the  details  only  for  the  area  bounded  by  the 
latitudes  i72°.  The  portions  enclosed  by  latitudes  72°4  <j>  4  84° 

-72°>  <f  >  -84° wi I  I  be  treated  summar i ly  later  on. 

At  stage  8  we  cut  the  octant  (truncated  at  y  -  *72°)  into  5  slices 
as  shown  in  figure  4.12(8).  The  meridians  are  eliminated  in  a  sequence 
implied  by  the  figure.  The  effort  results  from  the  following  table  4.4 


No . 

i 

j 

CPU 

1 

144 

684 

93.9 

2 

144 

456 

45.9 

3 

144 

532 

60.1 

4 

144 

456 

45.9 

Table  4.4.  Computational  effort 
for  the  4  phases  of  stage  8 

Summing  for  this  stage  we  obtain 


CPU  time  =  245.8  hours 


[stage  8] 


(4.11) 
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At  stage  7  we  have  40  blocks  in  the  shape  of  slices  as  shown  in  figure 
4.12(7).  Each  slice  is  split  into  4  blocks.  The  elimination  of  the 
dividing  lines  proceeds  in  a  sequence  as  shown  in  figure  4.12(7).  The 
effort  results  from  the  following  table  4.5 


No. 

i 

j 

CPU 

1 

30 

380 

26.6 

2 

30 

228 

10.  1 

3 

30 

228 

10.1 

Table  4.5.  Computational  effort 
for  the  3  phases  of  stage  7 

We  obtain  for  this  stage  a 

CPU  time  =  46.8  hours  [stage  11  (4.12) 

At  stage  6  we  are  dealing  with  160  18°xI8°  blocks.  We  split  one  block 
into  2  parts  by  a  spherical  surface  between  the  3rd  and  4th  layer. 
Confer  figure  4.12(6).  We  have  for  one  block 

i  =  4 

j  =  4*46  -  4*8  =  152 

The  resulting  CPU  for  all  160  blocks  is 

CPU  time  =  2.2  hours  [stage  6]  (4.13) 

We  ignore  the  upper  blocks  in  the  sequel.  Their  contribution  is 
marg i na I . 

At  stage  5  the  number  of  essent i a  I • b I ocks  is  still  160.  We 
decompose  anyone  of  them  into  3  subblocks  of  size  6  degrees  by  18 
degrees.  Cf .  figure  4.12(5).  There  are  two  phases,  and  each  of  them  has 

i  =  26 
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whereas 

j  =  4*38  -  4*5  +  2  =  134 
j  =  2*38  +  2*27  -  4*5  =  110 

The  effort  is 

CPU  time  =  21.8  hours  [stage  51  (4.14) 

At  stage  4  we  have  480  blocks.  (Figure  4.12(4)).  Each  one  is  subdivided 
into  three  6°x  6°  blocks.  Again  there  are  2  phases: 

i  =  6 

j  =  2*27  +  2*16  -  4*5  =  66 
j  =  2*38  +  2*16  -  4*5  =  88 

CPU  time  =  5.4  hours  [stage  4]  (4.15) 

At  stage  3  the  number  of  blocks  is  1440.  (Figure  4.12(3)).  Ue  eliminate 
the  4  nodes  at  the  lower  face  of  the  upper  subblock.  Consequently  we 
have  for  one  block: 

i  =  4 

j  =  4*16  -  4*5  =  44 

CPU  time  =  1.7  hours  [stage  31  (4  16) 

The  upper  portion  cut  off  in  this  way  does  not  have  any  interior  nodes 
and  is  therefore  ignored. 

At  stage  2  (see  figure  4.12(2))  we  eliminate  the  4  nodes  situated 
below  those  eliminated  in  the  previous  stage: 

i  =  4 

j  =  4*12  -  4*3  +  4  =  40 

.  CPU  time  =  1.5  hours  [stage  21  (4.17) 

At  the  first  stage  the  4  inner  nodes  at  the  bottom  of  each  of  the  1440 
blocks  are  eliminated  (see  figure  4.12(1)!) 


CPU  time  =  0.7  hours 


[s toge  i ] 


(4  I  8) 


Ue  now  return  to  the  chips  of  the  octants  cut  off  in  stage  9.  They  are 
located  between  ±72  and  ±84  degrees  of  latitudes.  Their  structure  is 
best  seen  from  figure  4.10(b).  Table  4.6  summarizes  the  elimination 
stages  necessary  to  decompose  these  blocks. 


Stage 

Dissecting  surface 

j 

j 

No  b 1 ocks 

CPU 

8' 

spherical  between 
3-rd  and  4-th  layer'" 

4 

304 

8 

0.4 

V 

A  =  const . , 

16 

264 

8 

1  .4 

producing  5  si  ices 

16 

132 

8 

0.4 

16 

176 

8 

0.6 

16 

132 

8 

0.4 

6' 

y  -  ±78° 

6 

81 

40 

0.2 

5' 

as  in  stage  3^ 

4 

44 

80 

0.  1 

4' 

as  in  stage  2 

4 

40 

80 

0.  1 

3' 

as  in  stage  1 

4 

28 

80 

0.0 

E 

-  3.6 

1)  Upper  portion  ignored  in  the  sequel 

2)  blocks  of  stage  5'  and  3  ore  identical.  They  are 
6°x6°  blocks  involving  layers  1,2,3. 


Table  4.6  Complementary  stages  for  regions  at  high  latitudes 


Hence  we  obtain 


100  - 


CPU  time  =  3.6  hours  [stages  3'  to  3'] 

The  total  CPU  is  obta  i ned  by  summing  over  all  stages 

Total  CPU  time  =  729  hours  (4.195 

Remark  •  We  have  been  somewhat  wasteful  by  a! lowing  many  elements  of 
small  width  (6  degrees  longitude)  at  latitudes  <f »  ±84  One  could  at  the 
lower  layers  imcgine  an  intermediate  ring  of  elements  between  84  and  86 
degrees  latitudes  whereby  the  longitudinal  width  wouid  be  18  degrees. 
This  would  probably  reduce  the  total  CPU  time  to  about  600  hours. 

In  any  case  it  turns  out  that  the  problem  is  not  managab 1 e  on  a 
computer  do 1 ng  only  half  a  million  steps  per  second.  However  there  are 
faster  computers.  The  ILLIAC  IV  is  described  in  Avila  et.al.  (1978)  as  a 
machine  that  has  64  processors,  each  of  a  speed  comparable  to  that  one 
assumed  above  All  processors  execute  the  same  instruction  at  a  time, 
but  each  one  operates  on  a  separate  data  stream.  Data  can  be  exchanged 
between  processors  (a  subset  of  the  processors  can  be  disabled  under 
program  control).  Such  a  machine  appears  to  be  well  suited  to  cut  down 
the  CPU  time  by  a  factor  which  may  approach  64.  The  elimination 
procedure  outlined  above  can  be  viewed  as  a  procedure  eliminating  nodes, 
whereby  each  node  contributes  8  parameters  to  the  system.  We  may  thus 
view  our  I  inear  system  as  one  having  as  many  unknowns  as  there  are 
nodes,  however  each  unknown  represents  actual ly  a  subvector  of  8 
elements.  From  the  viewpoint  of  eliminating  nodes,  an  elementary 
operation  is  thus  a  multiplication  of  two  8x8  matrices  followed  by  an 
addition.  There  should  be  a  way  to  organize  these  operations  such  that 
all  64  processors  of  the  ILLIAC  IV  are  busy  all  the  time. 

It  appears,  after  all,  that  a  global  solution  is  feasible  if  one 
of  the  worlds  best  computers  is  available.  However,  I  consider  it 
doubtful  that  such  a  computational  adventure  will  be  undertaken  in  the 
next  future.  Fortunately  enough  there  is  another  bonus  available  in 
physical  geodesy,  as  well  as  in  other  disciplines,  which  allows  us  to 
calculate  a  detai led  and  good  solution  in  a  local  area  of  interest 


without  being  forced  to  calculate  such  a  detailed  solution  everywhere. 
The  effect  which  we  are  going  to  exploit,  and  which  has  been  exploited 
a  lot  in  the  past,  is  called  "the  remote  zone  effect".  It  is  called,  by 
the  way,  "St.  Venant's"  principle  in  elasticity  theory. 

4.3.  The  remote  zone  effect. 

There  is  not  much  need  to  elaborate  on  it  at  length  Everybody  knows 
that  the  details  of  the  field  in  one  area  have  little  correlation  with 
the  details  in  remote  areas.  Put  it  in  other  words,  if  the  field  is 
changed  in  an  area  eg.  by  changing  the  mass  distribution  in  this  area, 
the  high  frequent  features  in  remote  areas  remain  practically  unchanged. 
In  applications  of  Stokes  and  Vening  Meinesz  formulas  this  effect  is 
utilized  in  a  way,  that  detailed  data  are  only  processed  in  a  fairly 
small  area  around  the  point  of  interest.  The  remote  zone  effect  can  also 
be  exploited  in  conjunction  with  other  methods,  eg.  the  surface  layer 
method,  the  collocation  method  and,  of  course,  the  finite  element 
method.  We  shall  consider  two  conf i gurat i ons .  The  first  one  is  a  strip 
in  which  a  detailed  geo  id  is  sought.  The  second  one  is  a  local  region, 
which  for  s i mp I i c i ty  w i I  I  be  taken  as  rectangularly  shaped. 

4.4,  Detailed  solution  in  a  strip 

Consider  a  partition  at  ground  level  as  shown  in  figure  4.13(a).  Assume 
that  this  pattern  extends  around  the  globe,  and  that  the  central  line 
corresponds  to  the  equator  The  size  of  the  smallest  squares  is  assumed 
to  be  2°x  2°.  The  partitioning  of  the  space  outside  the  surface 
(assumed  spherical)  is  demonstrated  in  fig.  4.13(b).  The  size  of  the 
elements  increases  as  we  go  away  from  the  surface.  The  outer  zone  is 
assumed  to  be  partitioned  into  specially  designed  elements  which  share 
two  more  nodes  at  the  north-  and  south-pole  See  also  fig.  4.14. 


Figure  4.13(a).  Partition  at  ground  level  of 


Figure  4.13(b).  Partition  of  vertical  cross 
section  AA'  as  indicated  in  figure  4.12(a). 
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Figure  4.14.  Element  partition  for  solution 
in  equatorial  strip 

The  potential  at  the  tuo  nodes  at  north-  and  south-pole  will  be  assumed 
to  be  known  Only  a  very  much  smoothed  version  of  the  potential  is 
needed  there,  Hence  these  nodes  will  not  contribute  to  the 
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computational  effort 

The  idea  is  to  get  a  detailed  field  in  the  vicinity  of  the  nodes 
located  at  a  vertical  plane  passing  through  the  axis  of  the  strip. 

At  the  final  stage  (stage  10),  we  eliminate  all  nodes  located  at 
the  main  profiles.  There  are  20  such  profiles,  each  having  60  nodes.  Ue 
obtain  a  system  which  is  structured  as  shown  in  figure  4.15 


Figure  4.15.  Cyclic  block  banded  structure  of  normal 
equations  at  the  last  stage  of  the  strip  solution. 

Ue  are  dealing  with  a  cyclically  blockbanded  system.  There  are  n  =  20 
diagonal  blocks.  The  size  of  each  block,  denoted  m,  results  from  the  60 
nodes  at  one  profile,  i.e.  m  =  60*8  =  480.  The  effect  to  tr i angu I ar i ze 
such  a  system  amounts  approximately  to 

-f-nni3  (4  20) 

steps.  Under  our  usual  assumptions  the  resulting  CPU  time  is 
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CPU  lime  =  3.28  hours  [stage  10] 


(4.21) 


At  the  next  lower  stage,  which  is  stage  9,  we  are  dealing  with  20 
blocks.  We  cut  each  into  two  parts  by  a  spherical  surface  indicated  in 
figure  4.13(b)  by  a  horizontal  arrow  labeled  9.  We  have  n  =  20  blocks, 
each  one  having  i  =  14,  j  =  2*60  =  120.  The  resulting  effort  amounts  to 


n 


ij  (i+J) 


(4  22) 


steps.  The 


l 

i 

j 

j 


CPU  lime  =  0.64  hours  [stage  9]  (4.23) 

No  further  effort  is  required  to  deal  with  the  upper  part  cut  off  at 
stage  9.  There  are  no  interior  nodes  left  in  it. 

At  the  next  stage  we  cut  off  two  portions  at  the  outside  of  any 
block  by  considering  two  vertical  cones  (surfaces  of  g>  =  il2°  and 
indicated  by  vertical  arrows  labeled  8  in  figure  4.13(b)).  We  have 
n  =  20  (number  of  blocks),  i  =  12,  j  =  2*54  =  108.  We  get 

CPU-time  =  0.44  hours  [stage  8]  (4.24) 

From  now  on  (stage  7)  we  deal  with  n  =  20  blocks,  each  one  having  a 
length  of  18°  and  a  profile  as  shown  in  figure  4.16 


4  4  4 

4  5  4 

Figure  4.16.  Vertical  block-profile  at  stage  7. 
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We  cut  off  the  upper  level  by  a  spherical  surface  indicated  by  an  arrow 
labeled  7  in  figure  4.16.  We  have  i  =  6,  j  =  2*40  +  2*11  =  102.  Hence 

CPU-time  =  0.19  hours  [stage  7]  (4.25) 

At  stage  6  we  cut  each  of  the  20  blocks  into  3  parts  by  2  planes 
vertical  to  the  ground  level  and  the  axis.  There  are  two  steps,  each  one 
having  i  =  26 .  The  number  of  junction  nodes  is  j  =79  and  j  =  70 
respect i ve I y . 

CPU-time  =  1.14  hours  [stage  63  (4.26) 

At  stage  5  we  are  dealing  with  n  =  60  blocks.  We  bisect  by  a  vertical 
axial  plane.  See  the  vertical  arrow  marked  5  in  figure  4.16.  We  count 
i  =  6,  j  =  2*35  =  70.  The  number  of  blocks  is  n  =  60.  Hence 

CPU  time  =  0.27  hours  [stage  53  (4.26a) 

At  stage  4  we  remove  the  lateral  portion  of  each  of  the  120  blocks.  See 
the  arrows  labeled  4  in  figure  4.16.  We  find  i  =6,  j  =  2*20  +  6  =  46 

CPU  time  *  0.25  hours  [stage  43  (4.27) 

The  profile  of  a  stage  3  block  looks  as  shown  in  figure  4.17.  The  depth 
of  such  a  block  is  6  degrees  of  longitude.  Hence  the  horizontal  cross 
section  is  square  shaped.  We  are  still  dealing  with  120  blocks.  At  each 
of  the  subsequent  3  stages  we  cut  by  a  spherical  surface.  See  the  arrows 
labeled  3,2,1  in  figure  4.17. 


Figure  4.17.  Vertical  b lock-prof i le  at  stage  3 
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At  stage  3  we  have:  i  =  4,  j  = 
CPU- time  =  0.12  hours 
At  stage  2  we  have:  i  =  4,  j  = 
CPU-time  =  0.10  hours 
At  stage  I  we  have:  i  =  4,  j  = 
CPU-time  =  0.03  hours 
The  total  CPU  time  is 

Total  CPU-time  =6.46 


4*14  -  4*4  =  40 
[stage  3] 
4*12  -  4*3  =  36 

[stage  2] 
4*8  -  4*3  =  20 
[stage  I] 

hours 


(4  28) 


(4.29) 


(4.30) 


(4.31) 


4,5.  Rectangular  region. 

Consider  a  region  of  interest  covering  an  area  of  32°x  64°  (the 
conterminous  United  Slates  are  contained  in  such  a  region).  Augment  this 
region  by  adding  layers  of  succesively  larger  elements  to  account  for 
the  remote  zone  effect.  The  element  design  at  ground  level  is  seen  from 
figure  4.18.  Note  that  the  element  size  at  surface  level  is  now  l°x  1°. 

The  vertical  design  is  based  on  4  layers  of  elements.  The  lower 
two  layers  follow  the  pattern  of  figure  4.18.  At  the  third  layer  we  fuse 
4  blocks  into  one,  and  we  do  the  same  at  the  fourth  layer.  The  basic 
building  block  of  our  element  design  is  thus  a  block,  which  we  call 
5-slage  block.  The  element  partition  at  one  of  its  four  vertical  faces 
is  shown  in  figure  4.19  (a).  There  are  25  nodes  on  a  face  and  7  nodes  at 
a  vertical  edge.  The  corresponding  blocks  for  the  augmented  zones  are 
shown  in  figure  4.19  (b),  (c) . 
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At  the  final  stage  (stage  12)  we  imagine  a  bisecting  vertical  plane 
passing  through  the  central  meridian.  See  the  arrow  labeled  12  in  figure 
4.18  We  count 

i  =  8*25  -  7*7  +  2*13  -  2*5  +  2*6  -  2*3  =  173 


Hence 

CPU  time  =  0.25  hours  [stage  12]  (4.32) 

At  stage  II  we  split  along  the  central  parallel.  Accordingly 
i  =  8*25  -  8*7  +  1*13  -  5  +  1*6  -  3  =  155 


whereas 


j  =  173 

There  a  two  blocks 

CPU  time  =  2.85  hours  [stage  II]  (4.33) 

At  stage  10  we  have  4  blocks.  Any  of  them  is  split  into  two  subblocks 

along  its  central  meridian; 

i  =  4*25  -  4*7  +  1*13  -  1*5  +  1*6  -  1*3  =  83 
j  =  12*25  -  11*7  +  2*8  +  2*3  =  245 

CPU  time  =  3.90  hours  [stage  10]  (4.34) 

At  stage  9  we  have  8  blocks.  We  split  by  parallels.  There  are  2  types  of 

blocks.  For  four  blocks  we  have 

i  =  4*25  -  5*7  =  65 
j  =  12*25  -  11*7  +  2*8  +  2*3  =  245 


CPU  lime  =  2.86  hours 


For  four  blocks  we  have 


i  =  4*25  -  4*7  +  1*8  +  1*3  =  83 
j  =  8*25  -  7*7  +  2*8  +  2*3  =  173 

CPU  time  =  2.20  hours 

Summing  for  this  stage  we  get 

CPU  time  =  5.06  hours  [stage  9]  (4.35) 

At  stage  8  we  split  again  by  meridians.  See  table  4.7 


No.  of  blocks  i  j  CPU  (hours) 


4 

29 

216 

0.88 

4 

47 

173 

1 .04 

4 

29 

209 

0.83 

4 

47 

137 

0.69 

Z 

=  3.44  hours 

Table  4.7. 

Contr i but 

ion  of 

var i ous 

block  types  at  stage  8 

We  see 

CPU  time  =  3.44  hours  [stage  8]  (4.36) 

In  the  following  stages  we  ignore  the  boundary  layers.  Our  estimation 
procedure  will  give  an  upper  bound  on  the  CPU. 

At  stage  7  we  assume  32  blocks  of  size  8° x  8°.  Splitting  iA  half 
means  i  =  29,  j  =  144 

CPU  time  =  3.32  hours  [stage  7]  (4.37) 

At  stage  6  ue  have  64  blocks,  i  =  II,  j  =  108 


1 1 1 


CPU  time  =  1.29  hours  [stage  6]  (4.38) 

At  stage  5  ue  have  128  blocks  of  size  4°x  4°.  We  remove  the  upper  layer: 
i  -  1,  j  =  72 

CPU  time  =  0.10  hours  [stage  5]  (4.39) 

At  stage  4  we  have  still  128  blocks.  We  bisect  by  a  vertical  plane:  i  = 
10,  j  =  4*21  -  4*5  +  1  =  65 

CPU  time  =  0.89  hours  [stage  4]  (4.40) 

At  stage  3  we  bisect  the  256  blocks  once  more:  i  =  3, 

j  =  2*16  +  2*8  =  48 

CPU  time  =  0.27  hours  [stage  3]  (4.41) 

Now  we  have  512  blocks  of  size  2° x  2.°.  We  remove  the  upper  half:  j  =  1 t 
j  =  4*8  =  32 

CPU  time  =  0.08  hours  [stage  22  (4.42) 

At  stage  1  we  remove  the  two  lowest  nodes  at  the  central  axis  of  a 

block:  i  =  2,  j  =  4*9-4*3  +  1  =  25 

CPU  time  =  0.10  hours  [stage  1]  (4.43) 

Summing  up  we  obtain; 

Total  CPU  time  =  21.55  hours  (4.44) 

This  is  still  quite  large.  On  the  other  hand,  if  one  owns  a  computer, 
such  computation  times  are  not  unrealistic. 


4,6.  Further  effort  to  cut  down  the  computation  time. 

Examining  the  above  models,  in  particular  the  global  one  and  the 
rectangular  one,  one  realizes  that  most  of  the  computation  time  is  spent 


I 


-112- 

at  the  higher  stages,  not  counting  the  final  stage.  The  large 
computation  time  is  partly  due  to  the  large  number  of  junction  nodes 
encountered  at  these  critical  stages. 

There  appears  to  be  a  possibility  to  further  cut  down  on  the 
computational  effort  by  more  sophisticated  use  of  the  remote  zone  effect 
(St.  Venants  principle)  Let  us  outline  the  procedure  on  hand  of  our 
last  example. 

Suppose  that  we  first  attempt  a  strip  solution,  but  one  with  an 
insufficient  number  of  elements  at  the  lateral  and  upper  part  of  the 
strip.  To  be  specific,  suppose  that  we  deal  with  the  strip  along  the 
central  parallel  in  figure  4.18  and  that  we  use  an  element  partition 
shown  in  figure  4.20(a)  (horizontal  partition  at  ground  level)  and 
figure  4.20(b)  (vertical  profile). 


(a)  (b) 

Figure  4.20 

Element  partition  for  auxiliary  strip  solution 
(a)  horizontal  section,  (b)  vertical  section 


At  the  final  stage  (stage  7)  we  deal  with  a  symmetric  b I ock-tr i d i agona I 
system.  In  contrast  to  figure  4.15  it  is  not  cyclical.  The  effort  for  n 
blocks  of  size  m  is  j-nm3.  In  our  case  n  =  1 6+ 1  =  17,  m  =  30*8.  (We  have 
increased  n  by  I  in  order  to  approx i mate  I y  account  for  the  contribution 
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of  the  boundary  elements).  Hence 

CPU  time  =  0.15  hours  [stage  7] 


(4.45) 


At  stcge  6  ue  cut  off  the  upper  blocks;  i  =  5,  j  =  60,  number  of  blocks 
=  16 


CPU  time  -  0.04  hours  [stage  63 


(4.46) 


At  stage  5  we  cut  into  halfes  by  a  vertical  plane  orthogonal  to  the 
axis;  i  =  22,  j  =  2*27  +  5  =  59,  number  of  blocks  =  16 


CPU  time  =  0.25  hours  [stage  5] 


(4.47) 


At  stage  4  we  remove  the  upper  blocks;  i  =  5,  j  =  54,  number  of  blocks 

32 


CPU  time  =  0.07  hours  [stage  43 


(4.48) 


At  stag'e  3  we  remove  the  lateral  elements.  We  have  i  =  4,  j  *  2*21  +  5 
47,  number  of  blocks  =  32 


CPU  time  =  0.04  hours  [stage  33 
We  have  32  blocks  shaped  as  in  figure  4.21 


(4.49) 


Figure  4.21.  Block  at  stage  2 

We  cut  into  halfes  by  o  vertical  plane  along  the  axis;  i  =  2,  j  =  39, 
number  of  blocks  =  32 


CPU  time  s  0.01  hours  [stage  23 


(4.50) 


At  stage  I  we  remove  the  two  remaining  inner  nodes;  i  =  2,  j  =  25, 


i 
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number  of  blocks  =  64 

CPU  time  =  0.0!  hours  [stage  I]  (4.5!) 

Hence 

Total  CPU-time  =  0.57  hours  (4.52) 

It  is  seen  that  such  a  profile  can  be  calculated  in  about  0.5  hours. 
Profiles  in  meridional  directon  require  about  0.25  hours. 


Figure  4.22.  Partition  of  rectangular  region  after 
auxiliary  profile  solutions  have  been  calculated 

Ue  calculate  3  long  and  7  small  profiles  as  shown  in  figure  4.22.  The 
CPU  time  =  3.72  hours 


The  potential  and  its  nodal  derivatives  are  not  thought  to  be  known 
precisely  at  the  axial  plane  after  a  profile  calculation  However  the 
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high  frequent  portion  is  thought  to  be  known  with  good  accuracy.  This 
means  that  in  the  further  treatment  of  the  rectangular  region  we  may  go 
along  with  a  smg I  I er  number  of  junction  nodes  at  the  vertical  planes 
along  profiles.  The  number  of  junction  nodes  along  a  profile  follows 
from  figure  4,19  Cb).  Considering  this  reduced  number,  the  computation 
time  for  the  various  stages  of  section  4.5  is  now  reduced  as  shown  in 
table  4.8. 


Stage  i 

J 

No . 

CPU 

12 

(final)  91 

1 

0  04 

II 

75 

91 

2 

0.36 

10 

43 

123 

4 

0.51 

9 

27 

123 

4 

0.29 

43 

91 

4 

0.31 

8 

II 

96 

4 

0.06 

27 

91 

4 

0.17 

11 

107 

4 

0  08 

27 

75 

4 

0.  12 

7 

29 

64 

32 

0.82 

6 

1 1 

66 

64 

0.51 

5 

(est imoled 

from  earl 

ier  values) 

<  0  10 

4 

10 

50 

128 

0  55 

3 

(estimated  from  earl 

ier  values) 

<  0.27 

2 

■ 

<  0.08 

1 

■ 

<  0.10 

E  <  4 .37 

Table  4.8.  Reduced  computational  effort  for  the  solution  in  a 
rectangular  region  after  the  calculation  of  auxiliary  profiles 
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We  obtain  a 

Total  CPU  time  =  8.09  hours 

Facit-  there  is  always  a  way  to  improve  on  a  previous  estimation.  We 
feel  that  also  the  700  hours  for  the  global  model  could  be  cut  by  the 
above  procedure  to  a  fraction  of  1/3  or  belter. 


5.  A  orooosal  for  a  numerical  solution  of  the  free  boundar 


value  problem 

5.1.  Isoparametric  elements. 

Isopar ame tr i c  elements  are  well  known  in  finite  element  analysis.  Confer 
e.g.  Strang-Fix  (1973),  section  3.3.  We  outline  the  basic  idea  on  hand 
of  a  two  dimensional  example.  Consider  the  r,  y  plane  and  consider  an 
element  partition  as  shown  in  figure  5.1(a).  Denote  by  SLJ(r,9>)  the  shape 
functons,  where  i  refers  to  the  node  and  j  to  a  specific  parameter  at 
this  node  (Cf.  section  3.2).  We  assume  C1  continuity  of  the  Stj(r,cp). 


F i gure  5 . I 

Mapping  from  straight  to  curved  elements. 

Consider  now  a  mapping  from  the  r,  y  plane  of  figure  5.1(a)  into  the 
R,  £  plane  of  figure  5.1(b).  The  mapping  functions  are 


R  «  R  ( r,  y) 
<J>  *  4>  ( r,  y) 


(5.1) 
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I 


fc 


>' 


The  elements  get  distorted  thereby.  The  distortions  sholl  be  locally 
smooth.  The  mapping  functions  shall  be  C 1  continuous.  In  fact,  we 
represent  them  in  terms  of  finitely  many  parameters  by  using  exactly  the 
shape  functions  StJ(r,<j»): 


R‘  ZRyS  (r.tf) 

ij  J  J 

(5.2) 

<P  * 

ij  J  J 


It  is  now  the  R,  $  plane  of  figure  5.1(b)  where  our  problem  in  terms  of  a 
differential  equation  of  the  field  and  in  terms  of  observations  on  this 
field  is  posed.  The  r,y>  plane  of  figure  5.1(a)  serves  only  an  auxiliary 
purpose  during  the  calculations. 

There  is  Laplace's  equation  for  the  field: 


3aV 

A V(R,<P)  • 


i  j_  TV  t 

R  <?R  dQ'1’ 


O 


(5.3) 


And  there  are  point  measurements  involving  V(R,<f>)  and  its  derivatives. 

In  order  to  have  something  specific  in  mind,  assume  that  the  image 
of  the  unit  circle  r  =  1  in  the  r,  y  plane  is  the  (finite  element 
approx i mat i on  of  the)  level  surface  of  the  reference  potential  in  the 
R, d  plane.  Assume  that  point  gravity  anomalies  are  prescribed  at  certain 
locat  ions  k  =  1,2,...,  K  which  are  the  images  of  r  *  I ,  .  Let  v- 

-  be  the  normal  to  the  level  surface  in  the  R,  4>  plane.  Then  the 
fo I  lowing  equat i ons  hold  (V  ...  disturbing  potential,  f  ...  normal 
grav i ty) : 


A9« 


±itv 

r  9*  R-RH.yJ 
♦  •  Qii,  vH) 


(5.4) 
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Our  next  problem  is  to  represent  the  potential  in  the  R,<i>  plane  in  terms 
of  finitely  many  parameters.  Here,  a  deliberate  detour  is  taken  One 
starts  representing  the  image  of  V  in  the  r,y  plane  in  terms  of  the 
shape  functions: 


V(r  y)  *  51  V  SLi(r.<f) 

tj  J  J 


(5.5) 


The  representation  of  V(R,4>)  is  then  obtained  by  means  of  the  inverse 


funct i ons 


r  •  r l R.<P)  •  R'1  (R&) 
s  ^  ( R .0)  1  <P~1(R,<P) 


(5.6) 


Fortunately  enough,  the  inverse  functions  will  never  be  used  explicitely 
during  the  numerical  calculations.  We  obtain 


V( R,<P)  *  Z  C  (r(R<*>)( 

ij  lJ  LJ 


(5.7) 


Remark :  The  fact  that  the  same  shape  functions  SLj  are  used  in  order  to 
represent  the  transformation  and  the  potential  is  responsible  for  the 
name  " i soparametr i c  elements". 

We  nou  outline  the  calculation  of  the  field  contribution  to  the 
normals.  Inserting  the  above  representat i on  of  V(R,d)  into  the  Laplacean 
(5.3),  ue  obtain 


Am, 


■•*)  *  Z-VrA  JR1  R  9R  Rl  9 <0*  J 


(5  8) 


The  following  manipulations  serve  to  circumvent  the  explicit  use  of 
r(R,<t>)/  gKR,4>)  during  the  subsequent  calculations.  Let  us  take  the  term 


FSij 

9Rl 


(5.9) 
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as  an  example.  It  holds  that 


3  Sa  , 

3Scj  3  r 

3^ 

.2*. 

3R 

3  r  9R 

9R 

3*5* 

9^  [  3  r 

3'Stj 

3RX 

3  ?7M 

/  + 

L  <?r3y 

3r 

9R 


3y 

30 


„  2Ll  <£l.  ,  i±±  d 

<9r  '  <9^ 


The  functions 

3^  3^;  3l^i  . 

ir; c5.ii) 

dont  require  any  further  treatment  because  the  SLj  are  simple 
expressions  in  r,  cf> ,  e  g.  bicubic  polynomials  uhich  are  differentiated 
as  shown  in  section  3.1.  Of  course,  after  differentiation  the 
substitution  r  =  r(R,40,  y  =  yCR,^)  has  to  be  imagined  but  this 
substitution  will  readily  be  undone  during  a  subsequent  change  of 
variables  in  the  integrals.  Only  the  derivatives  of  r,  <f  with  respect  to 
R  cause  trouble.  One  forms  the  Jacobian  matrix  of  partials 


3R 

3R 

3r 

3  y 

30 

3  r 

3>y 

(5. 12) 
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One  inverts  this  matrix' 


dr 

dr 

d<t> 

35 L 

djL 

_<?R 

3><p  _ 

1  i  n  th 

i  s  way 

dr 

is., 

SR  ' 

9R  * 

(5.13) 


(5.14) 


whereby  the  substitution  r  =  r(R,<$>),  y=  <f(R,<f>)  must  be  imagined,  but  is 
not  done  explicitely.  Denote  temporarily 
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Hence  the  second  derivatives  needed  in 


(5.18) 


are  calculated  In  an  analogous  fashion  all  quantities  needed  in 

eva  1  uat  i  ng  .dVCR, 4>)  can  be  obtained.  One  arrives  at  the  following  genera 

express i on 


&V(R,< P)  *  Z  VtjAS^trCRV  cfCw))  -- 

VciQ..ir(R,<t>)^(R<p)) 
ij  J  y 


(5.19) 


where  the  substitution  r  =  r(R/1^)^  <f  =  <f(R,4)  is  not  actually  carried 
out . 

We  proceed  to  the  integration  over  the  curved  element  in  the 
R, d  p 1  one : 


E  •t/' (AV(R.<t>))ZRcLRcL<t> 


RdRcL4> 


(5.20) 


Taking  the  variation  w  i  tTi  respect  to  the  V-.j.  one  gets  (as  in  section 
3.2.2.)  the  contribution  from  this  quad  to  the  normols  as 


£9, 


«-J 


v 


Kr  ■ 


(5.21) 


w  i  th 


Qc;  ( r(R.P),  <?(*,*))•  Q.,;,  ( rCR.ty'ft.R'*))  R dRd<t> 

u 


(5.22) 
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Now  comes  the  announced  transf ormat i on  of  parameters  of  the  r,  plane. 
The  domain  is  mapped  onto  the  uncurved  rectangle  in  the  r,  cp  plane 
Det[J(r,g>)]  is  the  Jacobian  determinant  of  this  transformation. 

giJ;Cy  :fQcj(r.y)  Qtf(r.tf)  R(r,<f)  ■  £et[j  (r.y)]  drdy>  (5  23) 

□ 

Remark  ■  Although  the  functions  QL](r  R(r,y)  readily  split  into  sums 
of  products  of  terms  F(r)G(cp),  the  Jacobian  Det[J(r,y)]  does  not.  Hence 
the  integration  procedure  is  not  easily  carried  out.  Numerical 
integration  may  be  necessary. 


Let  us  also  elaborate  on  the  contribution  of  the  measurements,  in  our 
special  case,  the  gravity  anomalies  First  we  need  a  representation  of 
the  vector  vh,  i  ,e.  the  normal  to  the  level  surface 


R  *  R(i,y) 
<p  *  <P  H  y) 


(5.24) 


We  may  view  this  as  a  parameter  r epresentat i on  of  the  surface.  In  the 
orthonormal  system 


e 


_  1  9 

e4>  ’  R  ' 


(5.25) 


the  vector  tangent  to  the  level 


surface  (curve)  has  the  r epresentat i on 


9fL 

dy 


(5.26) 


i 

i 


Hence  the  orthogonal  vector  has  the  repr esentat i on 


(5.27) 


(5.28) 


Discussion:  Remember  that  we  are  presently  dealing  with  a  completely 
known  isoparametric  transformat i on .  The  parameters  Ry,  ore 
prespec i f i ed .  What  is  the  benefit  from  such  a  transformation  and  what  is 
the  price  to  be  paid? 
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An  isoparametric  transformat i on  opens  the  way  to  exotica! ly  shaped 
reference  surfaces.  The  resulting  normal  equations  are  still  as  sparse 
as  those  for  the  spherical  Cor  ellipsoidal)  reference  surface  Hence  the 
solution  of  the  normals  require  no  additional  effort.  A  pr i ce  has  to  be 
paid  when  the  normals  are  formed.  The  integration  procedure  is  more 
costly  with  respect  to  programming  effort  and  computation  time.  On  the 
other  hand,  the  solution  of  the  normals  is  asymptotically  more  time 
consuming  than  the  formation.  If  N  is  the  No.  of  elements  at  ground 
level,  then  the  formation  of  normals  is  an  effort  of  cFN  where  on  the 
solution  is  an  effort  crN- .  Here  cF  and  cs  are  Cnearly)  constants. 
Recall,  by  the  way,  that  a  solution  by  means  of  surface  layer  elements 
or  collocation  requires  cN3  1  It  is  seen  that  i soparametr i c  elements 
affect  the  constant  factor  c,  but  not  cs  . 


5.2.  Approaching  the  free  boundary  value  problem  of  physical  geodesy 

Isoparametric  elements  as  outlined  in  the  previous  section  open  a  door 
to  approach  the  fundamental  problem  of  physical  geodesy  in  a  more  direct 
way  than  this  has  been  done  thus  far.  Our  subsequent  presentation  will 
be  expository.  Detailed  formulas  will  be  given  elsewhere. 

The  fundamental  problem  of  physical  geodesy  is  the  simultaneous 
determination  of  the  earths  figure  and  potential  from  measurements  of 
the  3-d i mens i ona I  gravity  vector  at  the  earths  surface.  The  historical 
approach  to  this  problem  seemingly  has  been  done  under  the  motto:  'First 
linearize  everything  in  sight  and  then  think  about  formulating  a 
meaningful  problem".  As  already  proposed  in  Meissl  (1971)  we  prefer  to 
formulate  a  problem  and  then  to  linearize  it.  However,  we  assume  that 
the  unknown  surface  of  the  earth  is  smooth.  This  means  that  the  terrain 
has  been  smoothed  and  that  the  gravity  measurements  have  been  corrected 
accord i ng I y . 
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Assume  that  a  sphere  of  radius  ro  situated  in  r,  y,  A  space  is  the 
pre-image  of  the  earths  surface  under  a  certain  unknown  transformation 

R  --  R(r.y./l) 

<P  °  &  Cr,<f,A)  (5.3!) 

A  *  A  ( r,  </,  A  ) 

Assume  that  the  measurements  g,,  gx,  ...  of  the  3  dimensional  gravity 
vector  are  taken  at  locations  identified  by  their  coordinates  rK  = 
r0 >  W*  /  AK  in  r,  y,  A  space  The  vectors  gK  themselves  are  considered  to 
be  represented  in  a  rectangular  equatorial  system.  Assume  further,  that 
the  earths  rotation  is  known  so  that  the  rotational  part  of  the 
potential  can  be  taken  into  account  by  a  known  reference  potential  U, 
while  the  unknown  disturbing  potential  V  can  be  assumed  harmonic. 

The  transformat  i  on  from  r,  y> ,  A  space  to  R,  <j>  ,  A  space  is  now  set 
up  according  to  equations  (5.2)  in  the  previous  section  as 

*-Z*SjSU(r.*A) 

ij  J 

<*>  '  Z>*  (5.32) 

ij  J 

A  *  HA.-  Sijiny.A) 

The  Ry,  dy  ,  ALj  are  assumed  known  for  all  nodes  i  except  those  situated 
at  the  earths  surface.  We  define  the  index  set  L  in  a  way  that  i e L 
identifies  precisely  these  nodes.  Thus  Ry,  ,  Ay  are  assumed  unknown 
for  i  6  L.  The  remaining  coefficients  Ry,  ALJ,  i^L  are  known  and  can 
be  chosen  according  to  computational  convenience.  For  example  they  can 
establish  a  transition  from  spherical  surfaces  r  =  const  to  ellipsoidal 
surfaces  in  the  R,  <t> ,  A  space. 
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Next  the  contribution  of  the  Field  to  the  normal  equations  is 
considered.  One  forms  the  energy  integral 

£  8  f  [  4>.A)  ]a  Rcos  <p  cLR  dL4>  cL/\  (5.33) 

Si. 

It  is  now  extended  over  the  entire  exterior  space  in  R,  <t>,  A  space 
Unknown  quantities  in  this  integral  are  not  only  the  nodal  values  Vt:  but 
also  Ry,  Alj,  ie  L.  Hence  variation  of  E  must  be  performed  with 
respect  to  all  these  quantities  Thereby  the  usual  decomposition  of  the 
domain  into  the  individual  elements  is  employed,  and  a  transformat i on  of 
the  element  integrals  back  to  r,  Cf>,  A  space  can  be  done  just  as  outlined 
in  the  previous  section.  The  resulting  normal  equations  are  linear  with 
respect  to  VLj  but  non  1  i  near  with  respect  to  Rp,  Ac ;j,  i  6  L.  They  are 
of  the  following  form 

£  9cj;  Cf  K'f  ’  °  <5.30 

The  equations  are  still  sparse  if  considered  as  linear  equations  in  Vy. 
However  sparseness  also  extends  to  the  nonlinear  equations  in  the 
following  sense.  Any  nonzero  coefficient  g cj;i]'  is  related  to  one  or  more 
elements  coupling  the  nodes  i  and  i'.  Hence  g^-cj1  will  only  involve  such 
transf ormat  i on  parameters  R„f,  0Ke,  AKe,  ksL  which  belong  to  these  common 
elements.  In  other  words,  if  the  nonlinear  normal  equations  are 
linearized,  then  the  resulting  equations  will  have  the  same  sparsity 
pattern  as  the  equations  of  the  earlier  chapters.  The  only  difference  is 
that  all  nodes  ieL  will  have  additional  parameters  Ry,  <t>c j,  Ay-. 

The  contribution  to  the  normals  from  the  assumed  gravity 
measurements  does  not  pose  much  difficulty.  The  observation  equations 
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ere  of  the  following  type 


§;{ui*.P.A) 

,  q(xl  ,  r(*J 

A) 

+  V(  R,<P,A)J 

*  9„W  *  r  w 

(5  35) 

+  V(fK,<t>,A)J 

■  9«"J  ^ 

U  is  the  reference  potential,  V  the  disturbing  potential  The  equations 
are  again  I  i  necr  in  V^,  but  noni  inear  in  <Pij,  Ay .  i  e  |_  The  reference 
potential  U  also  gives  a  contribution  due  to  the  unknown  location  of  the 
measurements .  i  .  e.  due  to  the  unknowns  R;,.  <P,; ,  Au  i  e  L.  The  residuals 

V  LJ  rj  • 

r**1,  rKCy^,  r„(1^  are  weighted  and  nonlinear  normals  are  formed  They  will 
show  an  analogous  sparsity  pattern  after  linearization. 

The  normals  of  field  and  measurements  are  added  end  solved  by 
Newtons  met-od.  This  method  requires  precisely  that  linearization  which 
we  were  talking  about  above. 

Ue  shall  conclude  this  outline  of  an  intended  research  project 
w i th  the  fo 1  lowing 

Remark i  It  is  not  considered  that  the  use  of  three-d i mens i ona I  gravity 
measurements  is  very  meaningful  in  the  geodetic  boundary  value  problem 
unless  additional  information  on  horizontal  position  is  available.  Such 
information  comes  either  from  ground  control  networks  or  from  space 
observations.  Hence  one  should  attempt  to  either  incorporate  this 
information  (without  destroying  the  sparsity  pattern)  or  to  remove  the 
horizontal  aegrees  of  freedom  in  the  unknown  mapping  from  r,  Cf ,  A  space 
to  R,  <P ,  A  space . 
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6.  Computer  experiments  for  2~d t mens i ono I  problems 
61.  Purpose  ana  scope 

The  computer  experiments  to  be  described  in  this  chapter  were  designed 
to  find  an  answer  to  the  following  specific  questions 

1)  Should  the  finite  element’  approach  be  based  on  the  Ritz-,  or  the 
Trefftz-,  or  on  the  least  squares  principle.  As  we  have  pointed  out  in 
section  3.3.1.  the  least  squares  principle  was  found  to  perform  best  in 
the  presence  of  noisy  and  redundant  data. 

2)  Are  cubic  polynomials  sufficiently  accurate,  or  should  quintic  poly¬ 
nomials  be  used?  The  question  is  perhaps  posed  in  an  overly  simplified 
way.  One  can  always  account  for  the  lower  degree  of  the  cubics  by 
choosing  a  smaller  element  partition.  Hence  the  question  should  be  asked 
as  follows.  Does  it  pay  off  to  replace  the  choice  of  cubics  and  an 
adequate  element  partition  by  the  use  of  quintics  and  an  element 
partition  having  appropriately  larger  elements?  Quintics  look  attractive 
because  they  are  C1  continuous.  One  can  even  enforce  the  Laplacean  to 
vanish  at  the  nodes.  Nevertheless  it  was  found  that  the  use  of  cubics  is 
prefer ab I e . 

3)  Can  the  at tenuat i on-w i th-a I t i tude-effec t  of  the  potential  be 
exploited  in  a  way  that  the  size  of  the  elements  increases  with  altitude 
fast  enough  to  ensure  that  the  total  number  of  nodes  is  bounded  by  a 
constant  times  the  number  of  nodes  at  the  surface.  The  answer  is 

aff irmat i ve. 

4)  What  is  the  best  way  to  represent  the  field  in  the  remote  outer  space 
of  the  earth?  It  is  believed  that  specially  designed  elements  of 
inifinite  size  together  with  appr opr i ate  I y  chosen  local  shape  functions 
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(different  from  cubic  polynomials)  are  the  best  choice.  Confer  section 

3.6 

5)  What  is  the  best  ratio  between  weights  applied  to  the  field 
contribution  to  the  normals  and  those  applied  to  the  contribution  from 
the  geodetic  observat i ons? 

6)  Is  there  a  way  to  combine  in  one  calculation  surface  data  with 
satellite  derived  spherical  harmonics? 

Due  to  the  large  CPU-times  predicted  for  3~d i mens i ona I  calculations  it 
was  decided  to  conduct  the  experiments  in  2  dimensions.  Stoke's  problem 
for  the  unit  circle  was  solved  by  means  of  finite  elements,  and  for  a 
set  of  artificially  generated  data.  Our  version  of  Stoke's  problem  is 
formulated  as  follows.  Find  a  potential  V(r,y)  in  the  outer  space  of  the 
unit  circle  such  that 

(>f)  V { r,oo)  =  0(75)  ,  r  -*  CO  (6.!) 

(•2 )  AV  --  0  )  r  >  1  (6.2) 

(3)  v(r,y)+  ~  »  f(<?)  ;  r  =  1  (6 .3) 

It  is  seen  that  we  have  eliminated  the  logarithmic  part  of  the 
potential.  V  may  be  viewed  as  a  disturbing  potential.  The  reference 
potential  U  may  absorb  the  logarithmic  part.  The  function  f(<y)  must  be 
free  of  "circular  harmonics’  of  the  zero-th  and  first  degree.  CI.e.  its 
Fourier  series  must  start  with  terms  in  cos  2f,  sin  2(f). 

Remark •  Obviously  the  stated  problem  is  most  easily  solved  by  means  of 
circular  harmonics,  i.e.  by  Fourier — analysis.  Representing  fCy)  as 

co 

cos  ny  *  dnsin  ny>) 

nil 


(6.3a) 
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the  solution  is  obtained  as 

V(r  y)  ■-  z_  p;  (  an  cos  ny  <•  bn  5in  ncf) 

n«2 


(6.3b) 


w  i  th 


O-n 


n  -  1 


(6.3 c) 


The  above  continuous  version  (61-3)  of  Stoke's  problem  is  replaced  by  a 
discrete  one.  First,  the  continuous  data  f(<f)  are  replaced  by  f(yK)  for 
discrete  arguments  Secondly,  the  potential  V(r,y)  is  replaced  by  a 
finite  element  representation 


X>cj  (6.4) 


Confer  sections  3.2.  and  3.3.!,  in  particular  equation  (3.19).  The 
element  partition  is  shown  in  principle  in  figures  3.6(a)-(b)  of  section 
3.3.2.  In  all  experiments  conducted  thus  far,  the  arguments  yK  were 
assumed  equally  spaced. 

The  outer  zone  was  represented  in  two  ways,  namely  by  a  "circular 
harmonics"  representation  of  comparatively  low  degree  (cf.  section  3.7, 
method  labeled  (4)),  and  alternatively  by  elements  of  infinite  size  as 
outlined  in  section  3.6. 

The  data  fCcp*)  were  generated  from  an  "assumed  potential"  of  the 

form 


V^sCr.cj)*  E  7^  l  a-n  cos  ny>  *  bns'n  ncp] 

n»  a 


wnth 


O-n 


(6.5) 
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The  coefficients  ot„,  fin  were  randomly  generated  in  the  interval  C — 0  S , 
+0.5].  The  maxima!  degree  N^,,  and  the  damping  exponent  d  were  varied. 
The  choice  d  =  1  corresponds  to  Kaula's  rule  of  the  thumb  (Equal  degree 
variances  of  the  radial  derivative  for  a  wide  range  2  C  n  €  NMAjr ). 

Remcrk :  Note  that  no  measurement  noise  was  assumed  to  be  superimposed 
upon  the  data  calculated  from  the  assumed  potential. 

The  procedure  described  in  chapter  3  and  section  4.1  yielded  the 
potential  at  the  nodes  of  the  chosen  element  partition  together  with  the 
nodal  derivatives.  This  ‘calculated  potential*  was  compared  with  the 
assumed  potential  (6.5).  Statistics  of  the  deviations  were  calculated 
and  tabulated.  The  calculations  were  carried  out  on  the  CSU-Computer 
Amdahl  470  V/6-II.  In  a  few  cases,  a  post-analysis  of  the  OSU  results 
was  done  on  a  desk-top  computer  WANG  2200  VP.  The  potential  was  then 
interpolated  into  the  interior  of  the  finite  elements  in  order  to  verify 
that  the  approximation  was  also  good  there.  Also  a  few  orbits  of  passive 
masspoints  were  numerically  integrated  for  the  assumed  potential  and, 
alternatively,  for  the  calculated  potential.  In  both  cases  a  circular 
symmetric  reference  potential  of  the  log  r  type  was  superimposed. 


6.2.  Parameters  dist inguishing  the  experiments. 


In  this  section  we  give  a  detailed  description  of  the  input  parameters 
characterizing  one  particular  experiment. 

IDEGR  degree  of  polynomials  in  elements  of  finite  size-' 

3  ...  bicubic  polynomials 
5  ...  biquintic  polynomials 
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ISUINF  a  switch  distinguishing  the  representation  of  the  field  in 

the  remote  outer  zone  r  >  rour 

1  ...  circular  harmonics 

2  ...  elements  of  infinite  size 


NLAYER  number  of  successive  layers  of  elements  in  the  circular 

ring  1  <  r  <  rgijr 

ITYPECI),  I  =  I,  ...,  NLAYER  type  of  elements  in  layer  number  I 

1  ...  simple  quads 

2  ...  compound  quads 


NELEM(I)  number  of  elements  in  layer  number  I  (the  lowermost) 

I  ayer 

HFCT  a  factor  governing  the  thickness  of  the  various  layers 

and  also  responsible  for  the  size  of  r0ur.  Cf.  the  following 
remark . 


Remark :  From  the  input  parameters  NLAYER,  ITYPECI);  1*1,  ...,  NLAYER, 
NELEM(I),  HFCT  the  element  partition  was  calculated  by  the  following  set 
of  formulas^ 


NELEMCI)  =  NELEMCI-I ) /ITYPECI- I );  1=2,  ...,  NLAYER  (6.6) 

DELPHICI)  =  2*  Jr/NELEMCI);  1=1,  ...,  NLAYER  (6.7) 

RADIUS ( I )  =  I 

RADIUS(I)  =  RADIUSCI-I )#(!  +  DELPHI(I)*HFCT) 

1=2,  . . .,  NALYER+!  (6.8) 
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Thereby  we  have  denoted 

NELEM(I)  number  of  elements  in  layer  I 

DELPHI(I)  angular  width  of  elements  in  layer  I 

RADIUS(I)  inner  radius  of  layer  I 
outer  radius  of  layer  I-! 

I  =  I ,  .  .  . ,  NLAYER 

RADIUSCNLAYER+ I )  =  rour  outer  radius  of  last  I ayer  =  outer  radius 
of  the  circular  ring  partitioned  into  elements  of  finite 
size. 

Ue  now  continue  to  describe  the  input  parameters. 

NPCTC  largest  degree  of  circular  harmonics  in  assumed  potential 

(denoted  NM-X  in  equation  (65)) 

DAMP  damping  factor  in  the  assumed  potential  (denoted  d  in 

equation  (65)) 

NPINST  this  parameter  governs  the  number  of  locations  at  which 

data  (gravity  anomalies)  were  calculated  from  the 
assumed  potential.  Remember  that  the  locations  are  equally 
spaced  NPINST  is  the  number  of  local ;ons  in  an  interval 
of  size  DELPHKi),  i.e  in  a  boundary  segment  of  a 
lowermost  finite  element  at  r  =  I .  The  data  are 
arranged  there  as  shown  in  figure  6.1.  It  is  seen  that 
the  interval  ends  correspond i ng  to  element  boundaries 
are  halfways  situated  between  two  measurements  locations. 
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The  total  number  of  fictitious  measurements  is  thus 
obtained  as  NELEMC ! )*NPINST .  Recall  that  no  simulated 
measurement  noise  was  superimposed  upon  the  data 


Figure  6.1.  Arrangement  of  fictitious  measurements  in 
a  boundary  segment  of  the  unit  sphere 


STKUGT  weight  applied  to  the  contribution  of  the  measurements 

to  the  normals.  (The  weight  for  the  field  contribution 
was  assumed  with  a  value  equal  to  1) 

The  following  4  parameters  apply  only  to  the  case  ISWINF  =  I,  i.e,  to 
the  representat i on  of  the  field  in  r  >  r^  by  circular  harmonics. 

NHARM  highest  degree  of  circular  harmonics  in  the  representat i on 

of  the  field  in  the  outer  zone.  This  representat ion 
i s  thus  g i ven  by 

Vtr.9)  *  H  — -  [  An  cos  ny>  ♦  Bnsin  nv)  (6  9) 
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The  calculated  values  for  the  unknowns  A„,  Bn  should 
approximate  those  of  a„ ,  bn  in  the  representation  (65) 
of  the  assumed  potential.  The  deviations  are  due  to 

(1)  discretization  of  data 

(2)  discretization  of  the  field  by  the  finite  element 
repr esentat i on 

(3)  NHARM  <  NPOTC 

(4)  roundoff  error 

this  parameter  governs  the  number  of  locations  at  uhich 
the  finite  element  representation  of  the  field  in  the 
ring  was  collocated  with  the  circular  harmonics 
representation  of  equation  (6.9)  in  the  outer  zone  The 
locations  are  all  found  at  r  =  rouT.  They  are  equally 
spaced.  There  are  NELEM(NLAYER)  intervals  at  r  =  r0(/r. 
NPINHA  gives  the  number  of  locations  in  one  of  these 
intervals.  Their  distribution  is  similar  as  that  shown  in 
figure  6.1,  i.e.  the  locations  are  equally  spaced  and 
the  interval  boundaries  were  assumed  halfways  between 
the  two  adjacent  locations.  The  total  number  of 
points  of  collocation  is  thus  NELEM(NLAYER)*NPINHA .  This 
number  was  frequently  assumed  appreciably  larger  than 
2*NHARM  -  2,  the  number  of  coefficients  in  (6.9).  This 
means  that  we  are  working  with  a  redundant  set  of 
locations  at  which  consictency  of  the  two  potentials 
was  enforced.  Hence  remark  3  given  in  section  3.7  applies 
mutatis  mutandis  (i.e.  translated  to  the  method  labeled 
(4)  in  section  3.7).  Consistency  can  only  be  required 
in  a  least  squares  sense.  Weight  assumptions  were 
necessary.  They  are  specified  by  the  subsequent  2 
weight  parameters. 


137  - 


HARWG0  weight  applied  to  the  values  of  the  potential  during 

collocation  at  the  locations  on  the  circle  r  =  r0ur 

HARWGI  weight  applied  to  the  radial  derivative  of  the  potential 

dur i ng  coll ocat ion.- 

Remark •  Recall  that  the  weight  applied  to  the  field  contribution  was 

assumed  equal  to  I . 

r 

The  following  3  parameters  apply  only  to  the  case  of  IDEGR  =5,  i.e.  to 

the  case  of  biquintics. 

N89  This  parameter  allows  to  choose  between  8 

parameters  per  node  and  9  parameters  per  node.  In 
case  of  8  parameters  per  node,  the  Laplacean  is 
fixed  to  zero  at  any  node.  This  gives  a  linear 
relation  between  the  nodal  parameters  Vrr  ,  Vr, 
which  was  used  in  order  to  eliminate  Vrr. 

FLPWGL  Because  it  was  observed  that  the  nodal  parameters 

involving  second  radial  derivatives  were  rather 
poorly  determined  at  r  =  I,  additional  fictitious 
observations  of  the  Laplacean  were  assumed  at  all 
measurement  locations.  The  appropriate  weight  was  FLPWGL. 

FLPUGU  This  input  parameter  is  analogous  to  the  previous 

one,  however  it  applies  to  r  =  rour. 
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6,3.  Detailed  results  for  two  experiments. 

6,3. i.  An  experiment  using  comparatively  large  elements. 


paramet< 

ers 

were 

specified  as  follows 

IDEGR 

= 

3 

[cub i cs! 

ISWINF 

= 

2 

[circular  harmonics  in  outer  zone] 

NLAYER 

= 

4 

[number  of  layers! 

ITYPECI) 

= 

1 

[simple  quads  at  first  layer] 

ITYPEC2) 

= 

1 

[simple  quads  at  second  layer] 

ITYPEC3) 

= 

2 

[compound  quads  at  third  layer] 

ITYPEC4) 

= 

2 

[compound  quads  at  fourth  layer! 

NELEMCt) 

= 

32 

[number  of  elements  in  first  layer! 

HFCT 

= 

1 

[approximately  square  shaped  elements! 

NPOTC 

= 

128 

[highest  degree  in  assumed  potential] 

DAMP 

= 

2 

[damping  factor  in  assumed  potential] 

NPINST 

= 

8 

[8  data  points  per  interval  at  r  =  1] 

STKWGT 

= 

5 

[weight  for  data] 

NHARM 

— 

8 

[highest  degree  of  circular  harmonics  in  the 
representat i on  of  the  field  in  the  outer  zone] 

NPINHA 

~ 

8 

[8  collocation  points  per  interval  at 

HARWG0 

- 

50 

[weight  in  collocating  V  at  r  -  rour  ] 

HARWGI 

= 

5 

[ueigth  in  collocating  Vr  at  r  =  rovrl 

Table  6.1  illustrates  the  geometry  of  the  element  partition.  A  pictorial 
representat ion  is  given  by  our  earlier  figures  3.6(a),  (b)  in  section 
3.3.2. 
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LAYER 

ITYPE 

NELEM 

RADIUS 

DELPHI 

1 

1 

32 

1 .00000 

0.19635 

2 

1 

32 

1  . 19635 

0. 19635 

3 

2 

16 

1  .43125 

0.39270 

4 

2 

8 

Hj  ur 

1  .99330 

3.55884 

0.78540 

Table  6.1.  Element  partition  in  experiment 
with  comparatively  large  elements 

The  assumed  potential  (equation  (6.5))  is  depicted  at  ground  level 
(r  =  I )  in  figure  6.2(a).  It  is  seen  that  the  potential  does  not  exceed 
the  value  0.14.  Recall  that  the  assumed  potential  was  used  to  generate 
the  discrete  set  of  gravity  anomalies.  From  these  data  the  potential  was 
calculated  backwards  by  the  finite  element  method.  The  calculated 
potential  at  r  3  1  is  shown  in  figure  6.2(b).  It  is  seen  that  the 
irregularities  are  somewhat  smoothed  out.  Figure  6.2(c)  shows  the 
difference.  Note  that  the  ordinates  are  now  scaled  differently.  We  see 
that  the  relative  error  is  about  2.5^. 

The  figure  6.3(a)-(c)  describe  in  a  similar  way  the  behaviour  of 
the  radial  derivative  Vr  at  r  =  I .  Figures  6.4(a)-(c)  are  devoted  to  the 
horizontal  derivative  Vy  at  r  =  I . 

It  is  obvious  that  the  accuracy  in  the  presently  described 
experiment  is  insufficient.  In  order  to  obtain  an  approx i mat  ion  of  the 
geo  id  at  the  cm  level,  the  finite  element  calculation  should  reproduce 
about  4  correct  digits  of  the  disturbing  potential  V.  Presently  only  2 
digits  are  correct.  We  have  nevertheless  exhibited  the  results  in  some 
detail  because  they  tell  us  very  instructively  what  we  can  expect 
qualitatively  from  a  finite  element  solution.  If  the  elements  are  chosen 
too  large,  the  details  of  the  field  can  not  be  properly  represented.  One 
should  nevertheless  expect  that  the  approximation  is  good  in  the  low 
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frequencies.  This  is  apparently  the  case  and  is  further  illustrated  by 
figures  6.2(d-e),  6.3(d-e),  6.5  and  6.6. 

Figure  6  2(d)  shows  the  assumed  field  V  at  r  =  I,  truncated  to 
circular  harmonics  of  degree  n  <  32.  This  low-frequent  part  of  the  field 
has  about  the  same  number  of  parameters  as  the  trace  of  the  finite 
element  representation  of  V  at  r  =  !.  The  two  graphs  6.2(b)  and  6.2(d) 
are  visually  nearly  indist inguishable.  Figure  6.2(e)  shows  the  high 
frequent  part  of  the  assumed  field,  composed  of  circular  harmonics  of 
degree  32  <  n  <  128.  The  two  graphs  of  figures  6.2(c)  and  6.2(e)  are 
aifferent,  but  the  magnitude  is  the  same.  Figures  6.3(d~e),  which  should 
be  compared  to  figures  6.3(b-c)  suggest  that  the  same  conclusion  holds 
for  the  radial  derivative;  the  finite  element  solution  is  about  as  good 
as  the  assumed  field  truncated  to  low-degree  harmonics  of  degree  n  <  32. 

Figure  6.5  shows  the  superposition  of  figures  6.3(a)  and  6.3(b). 
Figure  6.6  shows  the  calculated  Laplacean  AV  at  r  =  1 .  (The  assumed 
Laplacean  is  zero,  of  course).  The  calculated  Laplacean  is  discontinuous 
at  element  boundaries.  Hermite  bicubics  are  only  C1  continuous.  The 
graph  of  AV  suggests  that  the  finite  element  solution  achieves  a 
smoothing  by  shifting  -  in  a  balanced  way  -  positive  and  negative  masses 
outward  of  the  earths  body.  The  method  automatically  regularizes  the 
field  in  this  way.  The  shifting  of  masses  should  not  be  noti cable  at  the 
low  frequencies. 
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Figure  6.2(d).  The  assumed  field  V 
circular  harmonics  of  degree  n  <  3 
comDaral i ve I v  larae  elements. 


n 


a  I cu I  a  ted  minus  assumed 
i  in  experiment  using 


votive  VP  of  the  assumed  field 
cular  harmonics  of  degree  n  <  32. 
vely  large  elements. 


Figure  6.3(e).  Contribution  of  circular  harmonics  of 
degree  n  >  32  to  the  radial  derivative  Vr  of  the  assumed 
field  at  r  =  1.  Experiment  using  comparatively  large 
e I ements . 
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Figure  6.4(b).  Calculated  horizontal  derivative  Vy  at 
r  =  I  in  experiment  using  comparatively  large  elements. 
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In  order  to  put  the  potential  derived  from  finite  elements  to  a  further 
testy  a  number  of  orbits  of  passive  mass-points  was  numerically 
integrated.  A  reference  potential  equal  to 

U(r)  *  Log  -p  (6.10) 

was  assumed.  If  V^s^y)  and  VCLC(r,y)  denote  the  assumed  and  calculated 
potential  discussed  previously,  then  the  total  potential  used  in  the 
orbit  calculations  was  either 

K*SS(r.<j>Js  ^  +  (6.11a) 

or 

VCLC(r.y)  *  l/(r)  *  c0/ctc  (r.ij;  (6.11b) 

By  choos i ng 

Ca  *  0.0003 

an  attempt  was  made  to  relate  the  reference  potential  and  the  disturbing 


potential  in 

a  way  quant i 

tatively  similar  to  the  real 

earth. 

The  equation  of  mot 

ion  in  polar  coordinates  is 

•  m 

*  3 

dw 

r 

-  r g>4  » 

9r 

• 

jL 

(6.12) 

•  • 

n  r  y' 

4-  2  7*?  3 

r% 

For  near  circular  orbits 

one  may  represent  r  and  in 

the  fol lowing  way 

r » 

rQ  +  Ar 

(6.13) 
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Thereby  the  two  constants  r„,  cf0  fulfill 


ro% 


_  ^ 
-  i 


(6  14) 


One  proceeds  to  linearize  the  equation  (6  12)  obtaining: 


A r  -  %aAr  -  Zr0<j>0 Ay 
Ay  v  2~  Ar 

r  O 


-^T  *  Co-^tro+Ar^i'Acp) 

7— . •  4r  (ro+Ar;  %t  *  Ay) 

(r0  +  Ar)  % 


(6  15) 


We  shall  refer  to  C6.I2)  as  the  equation  of  the  ‘exact  orbit'  whereas 
(6.15)  are  the  equations  of  the  ‘differential  orbit". 

Either  version  of  the  orbital  equations  was  integrated  by  a  high 
degree  Runge-Kut ta-type  formula  found  in  Henrici  (1962),  p.  171.  When 
the  calculated  potential  was  used,  and  when  the  orbit  crossed  an  element 
boundary,  the  step  size  was  temporarily  decreased  by  a  factor  of  1/4. 

Figure  6.7(a)  shows  the  results  of  an  exact  orbit  calculation.  The 
orbit  passes  through  all  4  layers  of  the  circular  ring  subdivided  into 
elements  of  finite  size.  (Cf,  figure  3.6(b)).  The  radial  and  angular 
deviations  between  the  two  solutions  based  on  WCLC  and  W are  shown  in 
figure  6.7(b).  The  deviation  amounts  to  about  I  pp  300,000  This 
corresponds  to  a  20m-accuracy  in  the  real  world. 

Figures  6.8(a)-(b)  show  in  a  similar  way  the  results  for  a 
differential  orbit  with  respect  to  r0  =  I  .  I ,  =  0.90909...  It  is  a  low 

orbit  which  passes  about  midways  through  the  lowermost  layer.  We  see 
that  the  accuracy  is  about  I  pp  1,500,000  or  4-5  m. 

Remark  The  deviations  due  to  differing  Vct4  ans  V*fJ  of  the  orbits  (exact 
as  well  as  differential)  are  nearly  linear  in  c0  over  a  fairly  wide 
range.  Hence  the  assumption  c0  =  0.0003  can  be  changed  by  rescaling 
figures  6.7(b)  and  6.8(b)  appropr i ate  I y . 


Figure  6.7(a).  Exact  orbit  from  experiment  using 
compar a t i ve ! y  large  el emen t s .  Initial  va 1 ues  ar e 

p(0)  =  1.01,  F(0)  =  0,  y(0)  =  0,  9><0>  =  >  55. 

A  step  size  of  Ai  =  0.05  seconds  was  used  for  the 
calculated  potential  and  4t  =  0.0125  for  the  assumed 
potential.  A  value  c,  =  0.0003  was  chosen.  (Cf . 
equal i ons  (6 . I  1)5. 


0009010 


-.0000015 

-.0000020 


+.00004 


i 


+ .  90003 


+  .  00002 


+  .00001 


'.00000 


-.00003 


.00002 


-.00001 


+.00000 


-.00001 


-.00002 


-.00003 


-.00004 


-.00005 


-.00006 


-.00007 


-.00008 


-.00009 


Figure  6.8(a).  Radial  increme 
(bottom)  obtained  by  integrat 
assumed  potential.  Initial  va 
(0)  =  0.90909. . .  A  step  size 
potential,  and  =  C.0I25  f 
0.0003  was  chosen. 
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Remark •  Note  that  orbits  in  a  central  force  field  implied  by  (6.10)  are 
not  closed.  In  case  of  the  differential  orbit,  the  angle  between 
successive  peri  centers  (successive  apocenters)  is  given  by  (cf.  Arnold 
(1978),  p.  37,  problem  2) 


(6  16) 


Translated  into  time  this  gives 


-f-  *  sec  (6.17) 

<?o 


This  is  in  agreement  with  figure  6.8(a). 

We  summarize  some  further  results  of  this  experiment  in  tables  6.2 
and  6.3.  Table  6.2  lists  the  magnitude  of  the  nodal  parameters  and  the 
magnitude  of  the  deviations  between  calculated  and  assumed  field. 

Maximal  and  rms  quantities  are  listed  per  level.  There  are  7  levels  of 
nodes.  A  layer  of  simple  quads  is  bounded  by  two  levels.  A  layer  of 
compound  quads  has  one  additional  level  of  nodes  ha  If ways  between  the 
two  bounding  levels.  Level  1  corresponds  to  r  =  I,  level  7  to 
r  =  four  -  3 . 55884 .... 
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Level 

V 

Vr 

v. 

max 

rms 

max 

rms 

max 

rms 

! 

0. 126944 

0.066857 

0.502951 

0 

1 98303 

0.365901 

0.201368 

0.002493 

0.001005 

0.098707 

0 

043968 

0.092146 

0 . 039320 

2 

0.078590 

0.041004 

0.178782 

0 

093193 

0.212148 

0.1 13458 

0.000258 

0.000128 

0.005493 

0 

002216 

0.014307 

0.005119 

3 

0.049725 

0.025910 

0 . 08'7696 

0 

046367 

0.130977 

0.067720 

0  0001 36 

0.000055 

0.000715 

0 

000291 

0.002135 

0.000799 

4 

0.031793 

0.016990 

0.045800 

0 

024350 

0.075744 

0.043491 

0.000084 

0.000034 

0.000101 

0 

000054 

0 . 000380 

0.000221 

5 

0.021928 

0.011866 

0.026655 

0 

014127 

0.050927 

0.029458 

0.000093 

0.000033 

0.000075 

0 

000038 

0.000256 

0.000126 

6 

0.009972 

0.005744 

0.008395 

0 

004658 

0.020399 

0.014113 

0 . 000045 

0 . 000026 

0.000037 

0 

000024 

0.000612 

0.000418 

7 

0.005620 

0.003345 

0.003601 

0 

002040 

0.011373 

0.007977 

0. 000016 

0 . 000010 

0.000040 

0 

000022 

0.000126 

0 . 000090 

Table  6.2 

Maximal  and  rms  values  for  assumed  nodal  parameters 
(top  entry  in  each  field)  and  for  deviations  from 
calculated  values  (bottom  entries)  at  various  levels 


of  nodes.  Experiment  using  comparatively  large  elements. 
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Table  6.3  compares  some  of  the  assumed  potential  coefficients  a„,  bn; 
n  =  2,  NHARM  =  8  with  A„,  B„  obtained  by  the  collocation  procedure 

at  r  =  r„ur. 


A* 

Br> 

n 

a-» 

bo 

Oiff. 

Diff. 

-0.016589 

0.051241 

2 

-0.016752 

0.051570 

0.003163 

-0 . 000329 

-0.054197 

0.041 178 

3 

-0.055401 

0.042326 

0.001204 

-0.001148 

-0.000030 

0.015842 

4 

0.000313 

0.019369 

-0.000343 

-0.003527 

(the  remaining  coefficients  are  not  listed 
because  the  errors  are  comparable  in  size 
to  the  coefficients  themselves) 

Table  6.3 

Comparison  of  assumed  (top  entries)  and  calculated 
(middle  entries)  harmonic  coefficients.  Experiment 
with  comparatively  large  elements. 
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6.3,2.  An  experiment  using  comparatively  small  elements. 
The  inpul  parameters  were  specified  as  follows 


IDEGR 

= 

3 

[same  as  before] 

ISWINF 

= 

2 

[same  as  before] 

NLAYER 

= 

4 

[same  as  before] 

ITYPECI) 

= 

1 

[same  as  before] 

ITYPEC2) 

= 

1 

[same  as  before] 

ITYPEC3) 

= 

2 

[same  as  before] 

ITYPEC4) 

= 

2 

[same  as  before] 

NELEMC 1 ) 

= 

128 

[increased  by  a  factor  of 

4] 

HFCT 

= 

1 

[same  as  before] 

NPOTC 

= 

64 

[decreased  by  factor  1/2] 

DAMP 

= 

1 

[less  damping  than  before] 

NPINST 

= 

8 

[same  as  before] 

STKUGT 

= 

5 

Csame  as  before] 

NHARM 

= 

32 

[increased  by  a  factor  of 

4] 

NPINHA 

= 

8 

[same  as  before] 

HARWG0 

= 

50 

[same  as  before] 

HARWGI 

= 

5 

[same  as  before] 

Table  6.4  and  figure  6.9  illustrate  the  geometry  of  the  element 
partition.  Note  that  rovr  results  in  a  smaller  value  than  before. 


f 


1 
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LAYER 

ITYPE 

NELEM 

RADIUS 

DELPHI 

1 

1 

128 

! .00000 

0 . 04909 

2 

1 

128 

1 .04909 

0.04909 

3 

2 

64 

1  10058 

0.09817 

4 

2 

32 

t  .20863 

0.19635 

1.44595 

Table  1 

3.4.  El ement  par t i 

t i on  in  exper i ment  t 

comparat i ve I y  small  elements. 

Ue  now  exhibit  without  much  text  the  counterparts  of  figures  6.2-6. 
6. 7~6. 8  and  tables  6.2  and  6.3.  These  are  the  figures  6.10-6.14  and 
tables  6. 5-6. 6. 


Figure  6.9.  Element  partition  for  experiment 
using  comparatively  small  elements. 
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Figure  6.10.  Assumed  potential  V  Ctcp)  at  r  =  I  in  experiment  using 
comparatively  small  elements.  The  calculated  V  is  graphically 
indistinguishable.  The  bottom  figure  shows  the  enlarged  difference 
calculated  minus  assumed  V. 


♦2 


t'f'-f  L 


Figure  6. It.  Assumed  radial  derivative  Vr  Clop)  al  r  =  I  in  experiment 
using  comparatively  small  elements.  The  calculated  Vr  is  graphically 
nearly  indistinguishable  The  bottom  figure  shows  the  enlarged 
difference  calculated  minus  assumed  V-  at  r  =  1 . 


Figure  6.13(a).  Exact  orbit  in  experiment  using 
comparatively  small  elements.  Initial  values  are 

r(0)  =  1.01,  r  (0)  =  0,  g>(0)  =  0,  y(0)  =  1.15. 

A  step  size  of  0.0125  sec  was  used.  A  value 
c„  =  0.00003  was  chosen.  (Cf.  equations  (6.11)). 
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Figure  6.13(b).  Difference  in  r  (lop)  and  y  (bottom) 
of  exact  orbils  obtained  from  calculated  and  assumed 
potential.  Experiment  using  comparatively  small 
e I ements . 
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Note  that  the  assumed  potentials  are  different  in  the  examples  described 
in  the  previous  and  in  this  subsection.  In  subsection  6.3.1.  ue  had 
NPOTC  =  128  coefficients  in  the  assumed  potential  and  a  damping  factor 
of  DAMP  =  2.  In  this  subsection  NPOTC  =  64  and  DAMP  =  !.  The  gradients 
of  the  present  potential  are  about  10  times  larger  at  r  =  !  than  in  the 
earlier  example.  For  this  reason,  the  constant  c„  of  equation  (6. II)  was 
chosen  by  a  factor  of  t/10  smaller  than  before.  However,  a  remark  given 
in  section  6.3.1.  carries  over,  implying  that  the  figures  6.13(b)  and 
6.14(b)  may  be  rescaled  in  proportion  to  any  change  of  c0 . 

Due  to  the  differences  in  the  assumed  potentials  the  differences 
in  the  results  are  not  only  due  to  the  change  of  the  element  partition. 
Hence  the  results  are  not  immediately  compared.  In  the  next  section  a 
number  of  experiments  will  be  described  summarily.  From  this  additional 
information  it  may  be  concluded,  that  the  transition  from  NELEMO,)  =  32 
to  a  value  of  128  improves  the  results  by  about  I  digit  at  r  =  I. 

Figures  6.13  and  6.14  showing  the  results  of  the  orbit 
calculations  demonstrate  an  accuracy  of  about  I  pp  5,000,000  in  case  of 
the  exact  orbit.  This  corresponds  to  I -2  m .  The  differential  orbit  is 
good  to  I  pp  300,000,000.  This  correspond  to  2  cm. 
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Leve  1 

V 

— 

V, 

max  rms 

max  rms 

max  rms 

1 

0.560348  0.214973 

5.990022  2.260387 

6  250663  2.255706 

0.002427  0.000882 

0.046030  0.016975 

0328472  0.153894 

2 

0.384302  0.171249 

2.389440  0.797211 

2  324042  0  834182 

0.000698  0.000250 

0.053510  0.019825 

0.011507  0.004792 

3 

0.297927  0.143964 

1.330956  0.502973 

1  334445  0553725 

0.000094  0.000026 

0.004152  0.001534 

0  005652  0  001970 

4 

0.243919  0.123493 

0.853614  0.368318 

0.912405  0.427124 

0.000053  0.000023 

0.001245  0.000492 

0.005634  0.001950 

5 

0.209708  0.107151 

0.606866  0.286809 

0.696430  0  348634 

0.000044  0 . 000019 

0.000588  0.000227 

0  001248  0  000473 

6 

0.156166  0.082076 

0.348728  0.184487 

0  458199  0.248361 

0.000091  0.000046 

0.000463  0.000222 

0.003435  0.001451 

7 

0.122152  0.064510 

0.240170  0.127450 

0.340255  0.187136 

0.000058  0.000031 

0.000102  0.000052 

0.000789  0.000376 

Table  6.5 

Maximal  and  rms  values  for  assumed  nodal  parameters 
(top  entry  in  each  field)  and  for  deviations  from 
calculated  values  (bottom  entries)  at  various  levels 
of  nodes.  Experiment  using  comparatively  small  elements. 
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n 

A„ 
a, 
Diff . 

Bn 

bn 

Diff. 

1 - 

n 

i 

An 

On 

Diff. 

- — — - — 

Bn 

bn 

Diff. 

-0  033504 

0.103139 

-0.005873 

-0  026646 

2 

-0  033505 

0  103140 

13 

-0  006296 

-0  027436 

0 . y00yy 1 

-0.000001 

0.000423 

0 . 000790 

-0.166199 

0.126972 

0  009536 

0 . 027354 

3 

-0.166203 

0.126977 

14 

0.010171 

0.027769 

0 . 000004 

-0 . 000005 

-0.000635 

-0.000414 

0.001247 

0.077468 

-0.031018 

0  025260 

4 

0.001250 

0.077477 

15 

-0.030838 

8  025864 

-0.000003 

-0.000009 

-0.000180 

-0  000604 

-0  035952 

-0.058543 

-0.013724 

-0.010341 

5 

-0  035958 

-0.058558 

16 

-0  013725 

-001 1 124 

0 . 000006 

0.000015 

0.000001 

0.800783 

0.023489 

0.072583 

-0  007358 

-8.008160 

6 

0.023507 

0  072625 

17 

-0.007685 

-0.007283 

-0.000017 

-0 . 000042 

0  000247 

-0  000877 

-0.025322 

-0.013856 

-0.019404 

-0.015308 

7 

-0.025337 

-0  013859 

18 

-0.020553 

-0.014436 

0. 00001 S 

0.000083 

0.001149 

0.000872 

0.031356 

-0.051497 

0.022968 

-0.018957 

8 

0.031431 

-0.051597 

19 

0.023800 

-0.020617 

-0.000075 

0.000100 

-0.000832 

0.001660 

0.052027 

-0.032417 

-0.012259 

-0  015923 

9 

0.052222 

-0.032523 

20 

-0.012340 

-0.016005 

-0.000195 

0.000106 

0.000081 

0.000082 

0.003373 

0.01737! 

-0.001372 

0.007838 

10 

0.003460 

0.017520 

21 

-0.001379 

0.007600 

-0.000087 

-0.000149 

0.000006 

0.000238 

0.012190 

0.018735. 

-0.020856 

0.013256 

1 1 

0.012312 

0.018955 

22 

-0.019737 

0.011703 

-0.000122 

-0.000220 

-0.001 119 

0.001553 

0.000715 

-0.001960 

-0.004177 

0.005744 

12 

0.000830 

-0.002129 

23 

0.000844 

0.608234 

-0.000016 

0.000169 

-0.005021 

-0.002490 

Table  6.6 

Comparison  of  assumed  (top  entries)  and  calculated 
(middle  entries)  harmonic  coefficients.  Experiment 
using  comparatively  small  elements. 
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5.4.  Summary  of  other  experiments. 

6.4.1.  Experiments  using  bicubics 

Table  6.7  lists  input  parameters  and  resulting  accuracies  for  a  number 
of  exper i men t s  using  b i cub i cs . 


Exper i ment 

1 

2 

3 

4 

5 

6 

7 

8 

9 

ISUINF 

1 

1 

1 

2 

1 

2 

1 

1 

1 

NLAYER 

4 

4 

4 

4 

4 

4 

5 

4 

4 

ITYPEC 1 ) 

1 

! 

I 

1 

1 

1 

1 

1 

1 

ITYPEC2) 

1 

1 

1 

1 

1 

1 

1 

1 

1 

ITYPEC3) 

2 

2 

2 

2 

2 

2 

2 

2 

2 

ITYPEC4) 

2 

2 

2 

2 

2 

2 

2 

2 

2 

ITYPEC5) 

/ 

/ 

/ 

/ 

/ 

/ 

2 

/ 

/ 

NELEMC!) 

32 

32 

32 

32 

32 

64 

64 

128 

128 

HFCT 

1 

1 

1 

1 

1 

1 

0.75 

1 

1 

rour 

3.56 

3.56 

3.56 

3.56 

3.56 

2.01 

2.72 

1.45 

1  .45 

NPOTC 

16 

32 

128 

128 

128 

32 

64. 

64 

128 

DAMP 

1 

1 

2 

2 

1 

1 

1 

1 

1 

NPINST 

8 

8 

8 

8 

8 

8 

8 

8 

8 

STKWGT 

5 

5 

5 

5 

5 

25 

5 

5 

5 

NHARM 

8 

8 

8 

/ 

8 

/ 

8 

32 

32 

NPINHA 

8 

8 

8 

/ 

8 

/ 

8 

8 

8 

HARWG0 

50 

50 

50 

/ 

50 

/ 

50 

50 

50 

HARUG! 

15 

15 

5 

/ 

5 

/ 

5 

5 

5 

(table  continued  on  next  page) 


Exper i ment 


IE-3  2E-3  3E-4 
4E-3  5E-3  IE-3 
IE-2  2E-2  4E-3 


/  2E-3 
/  6E-3 
/  2E-2 


8  9 


8E-4  IE-6  IE-5 
2E-3  5E-6  6E-6 
IE-2  9E-6  8E-6 


IE-3  2E-2 

IE-3 

IE-3 

4E-2 

IE-3 

IE-2 

9E-4 

7E-3  i 

5E-4 

3E-3 

IE-4 

IE-4 

3E-3 

4E-4 

8E-4 

2E-4 

IE-3  ! 

IE-4 

IE-3 

5E-5 

6E-5 

IE-3 

5E-5 

5E-4 

2E-5 

3E-4 

8E-5 

7E-4 

3E-5 

4E-5 

7E-4 

4E-5 

4E-4 

2E-5 

IE-4 

9E-5  5E-4 

3E-5 

4E-5 

5E-4 

4E-5 

3E-4 

2E-5  5E-5 

5E-5 

4E-4 

2E-5 

IE-5 

4E-4 

7E-5 

2E-4 

5E-5 

5E-5 

9E-6 

2E-4 

9E-6 

IE-5 

3E-4 

7E-5 

2E-4 

3E-5 

3E-5 

/ 

/ 

/ 

/ 

/ 

/ 

2E-4 

/ 

/ 

/ 

/ 

/ 

/ 

/ 

/ 

9E-5 

/ 

/ 

7E-3  5E-1  4E-2 

4E-2 

3 

IE-2 

7E-I 

2E-2 

1 

8E-3  5E-2  2E-3 

2E-3 

5E-2 

IE-2 

2E-2 

2E-2 

IE-1 

6E-4  6E-3  3E-4 

3E-4 

7E-3 

IE-3 

5E-3 

2E-3 

IE-2 

2E-4  IE-3  5E-5 

5E-5 

IE-3 

3E-4 

2E-3 

5E-4 

2E-3 

IE-4  SE-4  4E-5 

4E-5 

5E-4 

2E-4 

IE-3 

2E-4 

7E-4 

6E-5  2E-4  2E-5 

3E-5 

3E-4 

9E-5 

7E-4 

2E-4 

2E-4 

5E-5  2E-4  2E-5 

IE-5 

2E-4 

2E-4  5E-4 

5E-5 

7E-5 

/  /  / 

/ 

/ 

/ 

3E-4 

/ 

/ 

/  /  / 

/ 

/ 

/ 

2E-4 

/ 

/ 

6E-2 

8E-! 

4E-2 

4E-2 

2 

IE-1 

9E-I 

2E-I 

1 

9E-4 

IE-1 

5E-3 

5E-3 

2E-1 

3E-3 

6E-2 

5E-3 

2E-I 

2E-3 

2E-2 

8E-4 

8E-4 

2E-2 

2E-3 

IE-2 

2E-3 

2E-3 

IE-3  2E-3 

2E-4 

2E-4 

3E-3 

IE-3 

6E-3 

2E-3 

3E-3 

5E-4 

IE-3 

IE-4 

IE-4 

IE-3 

3E-4 

2E-3 

5E-4 

8E-4 

IE-3 

IE-3 

4E-4 

3E-4 

2E-3 

IE-3 

3E-3 

2E-3 

2E-3 

3E-4 

6E-4 

8E-5 

2E-4 

7E-4 

3E-4 

6E-4 

4E-4 

4E-4 

/ 

/ 

/ 

/ 

/ 

/ 

2E-3 

/ 

/ 

/ 

/ 

/ 

/ 

/ 

/ 

9E-4 

/ 

/ 

Table  6.7 

Summary  of  experiments  using  bicubics.  Experiment 
No.  3  and  No.  8  are  those  discussed  in  detail  in 
sections  6.3.1.  and  6.3.2. 
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6.4.2.  Experiments  using  biquintics. 

Table  6.8  lists  input  parameters  and  resulting  accuracies  for  a  number 
cf  experiments  using  biquintics.  The  parameters  NLAYER  (=4),  ITYPE(I) 
(=1),  ITYPEC2)  <=!),  ITYPEC3)  (=2),  ITYPEC4)  (=2)  are  not  repeatedly 
listed  because  the  indicated  values  were  the  same  for  all  experiments 
documented  in  this  table. 


Exper i ment 

1 

2 

3 

4 

5 

NELEMC1) 

32 

32 

32 

32 

32 

HFCT 

1 

0.75 

0.75 

1 

1 

p 

'  ovr 

3.56 

2.71 

2.71 

3.56 

3.56 

NPOTC 

32 

128 

128 

128 

128 

DAMP 

1 

2 

1 

1 

1 

NPINST 

8 

16 

16 

16 

8 

STKWGT 

1 

1 

1 

1 

1 

NHARM 

8 

8 

8 

8 

8 

NPINHA 

8 

8 

8 

8 

8 

HARUG0 

10 

10 

10 

10 

10 

HARUGI 

2.5 

2.5 

2.5 

2.5 

2.5 

N89 

9 

9 

9 

9 

8 

FLPWGL 

IE-3 

IE-3 

IE-3 

IE-3 

IE-3 

FLPWGU 

3E-3 

IE-4 

IE-4 

IE-4 

IE-4 

(table  continued  on  next  page) 


accuracy  af  levels 
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Exper i ment 

1 

2 

3 

4 

5 

2 

IE-3 

1E-S 

4E-4 

5E-4 

3E-2 

A,  B 

3 

2E-4 

4E-6 

2E-4 

3E-4 

2E-3 

4 

3E-4 

3E-5 

4E-4 

6E-4 

5E-3 

1 

2E-3 

4E-4 

3E-2 

3E-2 

3E-2 

2 

2E-3 

7E-6 

2E-4 

5E-4 

2E-3 

3 

IE-3 

7E-6 

2E-4 

3E-4 

IE-3 

V 

4 

9E-4 

5E-6 

2E-4 

2E-4 

7E-4 

5 

6E-4 

4E-6 

IE-4 

IE-4 

5E-4 

6 

2E-4 

2E-6 

6E-5 

5E-5 

2E-4 

7 

9E-5 

IE-6 

4E-5 

3E-5 

2E-3 

</» 

QJ 

> 

1 

IE-2 

3E-2 

2 

2 

2 

QJ 

2 

9E-3 

3E-4 

IE-2 

3E-2 

2E-2 

U 

3 

3E-3 

4E-5 

2E-3 

5E-3 

6E-3 

<_ 

a 

O 

Vr 

4 

IE-3 

8E-6 

3E-4 

2E-4 

9E-4 

03 

5 

9E-4 

6E-6 

2E-4 

2E-4 

6E-4 

6 

3E-4 

6E-6 

8E-5 

6E-5 

2E-4 

7 

6E-5 

2E-6 

3E-5 

6E-5 

IE-4 

1 

4E-2 

3E-2 

2 

2 

2 

2 

IE-3 

IE-4 

4E-3 

7E-3 

9E-3 

3 

2E-3 

2E-5 

9E-4 

IE-3 

2E-3 

vy 

4 

IE-3 

5E-5 

5E-4 

3E-4 

IE-3 

5 

8E-4 

2E-5 

3E-4 

2E-4 

9E-4 

6 

3E-4 

8E-5 

2E-4 

IE-4 

4E- 4 

7 

2E-4 

2E-5 

IE-4 

2E-5 

3E-3 

Table  68.  Summary  of  experiments  using  biquinlics. 
Add  it  ional  input  parameters  were  NLAYER  =  4, 
ITYPECI)  =  I,  ITYPEC2)  =  I ,  ITYPEC3)  =  2,  ITYPEC4) 
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6.4.3.  Pi scuss i on . 

II  is  rnt  eresting  to  compare  experiments  based  on  the  same  assumed 
fie! d,  e.g.  the  bicubic-experiments  5  and  9.  (The  field  at  r  =  !  is 
shown  in  figures  6.10-6.12.  Also  table  6.5  applies  to  this  field)  It  is 
seen  that  the  transition  from  32  elements  in  the  lowermost  layer  to  128 
elements  results  in  an  increase  of  accuracy  by  one  digit.  On  the  other 
hand,  experiments  8  and  9  differ  only  with  respect  to  NPOTC.  One 
recognizes  that  experiment  9  looses  accuracy  in  the  lower  layers  due  to 
a  rougher  field  there.  However  the  accuracy  in  the  upper  layers  is  about 
the  same.  This  confirms  our  reasoning  that  the  finite  element  solution 
regularizes  the  field  without  affecting  the  lower  frequencies.  A 
comparison  of  experiments  3  and  4  shows  that  the  representation  of  the 
field  in  the  outer  zone  by  spherical  harmonics  is  about  equivalent  to 
representation  by  elements  extending  to  infinity.  However  it  must  be 
borne  in  mind  that  the  computational  effort  is  less  in  the  case  of 
elements  extending  to  infinity. 

Experiment  No.  5  using  bicubics  and  experiment  No.  4  using 
biquintics  are  otherwise  based  mostly  on  equal  parameters.  A  comparison 
shows  that  the  improvement  coming  from  biquintics  is  marginal. 

Experiments  No.  4  and  No.  5  of  the  b i qu i nt i cs-tab I e  show  that 
fixing  the  Laplacean  at  zero  gave  worse  results.  It  is  particularly 
interesting  how  poorly  the  harmonic  coefficients  were  recovered  in 
experiment  5.  In  all  experiments  using  biquintics  the  recoverage  of  the 
harmonic  coefficients  was  inferior  to  that  in  the  experiments  using 
bicubics.  The  reason  is  not  yet  completely  clarified. 


The  programs  for  the  Amdahl  470  V/6— II  were  written  in  an  extended 
FORTRAN.  They  comprised  about  6000  statements.  Double  precision  was  used 
throughout.  Coding,  testing  and  processing  required  about  6  weeks.  No 
attempt  was  made  to  optimize  speed  by  using  assembly  code  in  the  inner 
loop  of  the  equation  solver.  The  CPU  time  for  the  largest  experiment  was 
83  seconds.  It  required  the  solution  of  a  linear  system  with  2368 
unknowns.  Recherches  show  that  about  half  of  the  time  was  spent  on 
evaluating  the  assumed  potential  and  its  derivatives  for  comparison 
purposes.  The  time  spent  on  the  solution  of  the  linear  system  must  have 
been  less  than  40  seconds.  A  full  system  of  this  size  would  have 
required  a  CPU  of  more  than  74  minutes.  Hence  a  factor  of  about  1/100 
has  been  gained  by  the  nested  dissection  method. 

The  programs  for  the  post-analysis  on  the  WANG  2200  VP-desk-top 
computer  were  written  in  an  extended  BASIC,  called  WANG  BASIC-2.  About 
500k  bytes  of  code  were  assembled  requiring  an  effort  of  about  3  weeks. 
Doing  the  post-analysis  on  a  large  computer  would  have  required  a  time 
span  of  probably  3  months.  Small  desk  top  computers  have  powerful 
editing  facilities,  they  offer  instant  response  during  editing,  and 
instant  diagnosis  during  testing.  It  is  a  nonsense  to  solve  small 
problems  on  large  computers. 


6.4.5.  Further  desirable  experiments. 

Time  pressure  did  not  allow  to  conduct  all  experiments  the  author  had 
originally  in  mind.  It  was  intended  to  generate  "assumed  potentials' 
other  than  those  defined  by  circular  harmonics.  For  example  it  would  be 
interesting  to  see  how  the  method  performs  for  a  potential  generated  by 
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buried  dipoles.  Also  potentials  with  other  types  of  near  singularities 
could  be  tried  Another  test  could  have  been  devoted  to  irregularly 
distributed  measurement  locations,  to  noisy  data,  and  to  less  regular 
element  partitions  Comparisons  with  the  results  of  other  methods 
applied  to  the  same  2-dimensional  problem  are  missing.  Finally,  tests  in 
3  dimensions  should  be  performed.  Unfortunately  there  is  a  limit  to  the 
cnount  of  work  a  single  person  can  do  in  a  year,  in  addition  to  teaching 
and  administrating.  It  is  intended  to  continue  the  experiments  with  the 
help  of  my  coworkers  and  students  and  to  give  a  more  complete  report  at 
a  later  occasion 
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7 ,  A  proposed  hybrid  method. 

7.1.  Revision  of  the  surface  layer  method 


There  is  an  obvious  way  to  introduce  finite  elements  to  the  surface 
1 ayer  method,  namely  to  represent  the  surface  layer  density  in  terms  of 
2~6 i mens i ona I  finite  elements.  Bilinear  functions  defined  over 
rectangles  in  the  if,X  plane  would  ensure  a  continuous  surface  layer 
density.  Bicubics  would  give  a  C'  function.  No  improvement  in 
computational  efficiency  can  be  expected  from  such  a  procedure.  The 
normal  equations  would  still  be  fully  occupied.  The  number  of  parameters 
per  node  increases  from  one  to  four.  Hence  block  areas  can  be  chosen  4 
times  as  large  as  in  the  case  of  constant  densities  within  blocks.  The 
benefit  would  be  a  field  having  continuous  derivatives  down  to  the 
surface  of  computation. 


7.2.  Mu  I t i po I e  I ayer  . 


Keeping  in  mind  that  the  potential  V  to  be  represented  is  actually  a 
disturbing  potential,  and  that  frequently  the  reference  potential  U  is 
chosen  in  a  way  that  V  is  free  of  harmonics  up  to  a  certain  degree  N-l, 
it  is  tempting  to  try  a  surface  layer  density  such  that  the  generated 
potential  is  free  of  harmonics  up  to  degree  N-l.  Such  densities  are 
available,  though  not  in  the  form  of  a  single  layer.  Layers  of 
multipoles  must  be  chosen.  The  potential  generated  by  a  dipole  is 

V  =  fLcosJ  (7. ,5 

l 

Here  fi  denotes  the  strength  of  the  dipole.  The  meaning  of  y  and  l  is 
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point  at  which 
V  is  evaluated 


Figure  7.1.  Explains  the  notation  of 
quantities  related  to  a  dipole 

It  is  useful  to  remember  that  the  potential  of  a  dipole  is  obtained  from 
that  one  of  a  mass  point  by  taking  the  directional  derivative  with 
respect  to  the  axis  of  the  dipole  and  at  the  location  of  it. 

cos  p  3_  I  (7  2) 

pr"  ~  dv  r 

It  follows  that  the  potential  of  a  dipole  decreases  like  as  r— ®. 

Hence  it  is  free  of  zero  and  first  order  harmonics.  A  dipole  layer  is 
obtained  by  locating  dipoles  at  a  surface  and  letting  the  strength  of 
the  dipole  be  a  function  of  location.  The  axis  may  also  vary  with 
location.  Physically  and  mathematically  most  meaningful  is  the 
coincidence  of  the  axis  with  the  surface  normal.  The  potential  of  a 
dipole  layer  is  also  0(^7} . 

Dipoles  may  be  generalized  to  multipoles.  Multi  poles  were  already 
studied  by  Gauss  and  Maxwell.  An  N+l  pole  of  unit  strength  is  obtained 
as 


seen  from  figure  7.1 


J _ <L  .  ..  .  JL  1  (73) 

Bv,  Bv2  3vn  l 

Confer  e.g.  Lense  (1954),  p.  80  ff.  Trying  most  simple  things  first,  one 
puts  V|=V2=  ...  =vN=v,  thus  letting  the  N  axes  coincide  and  obtaining 
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I 


i 


i 


3n  I 

z7  { 


(7.4) 


The  generated  potential  is  o[r~jp;  i.e.  it  is  free  of  spherical 
harmonics  up  to  and  including  degree  N~1 .  The  same  holds  for  the 
potential  of  a  multipole  layer.  The  axes  may  be  chosen  in  coincidence 
with  the  surface  normal. 

What  is  the  advantage  of  taking  a  mullipole  layer  instead  of  a 
single  layer?  At  a  first  glance  it  appears  that  mullipole  layer  also 
generates  a  full  normal  equation  matrix.  Mind,  however,  that  the  kernel 
of  a  mullipole  of  order  N+l  decreases  as  0{r  n  »'i }  as  r— ' This  implies 
that  even  for  moderately  large  N  the  contribution  of  the  multipoles 
located  at  a  small  surface  element  decreases  rapidly  as  one  moves  away 
from  that  surface  element.  If  the  surface  on  which  the  multipole  layer 
is  assumed  is  subdivided  into  finite  elements,  and  if  the  strength 
function  is  represented  in  terms  of  nodal  parameters,  then  the  nodal 
parameters  of  elements  at  a  larger  distance  will  be  coupled  only  by  very 
small  coefficients  which  may  be  replaced  by  zeroes.  Hence  a  sparse 
system  of  norma!  equations  will  be  obtained. 

What  are  the  difficulties  to  be  expected  in  implementing  the 
multipole  layer  method,  and  what  is  an  appropriate  finite  element 
representation  of  the  strength  function?  How  smooth  should  it  be  chosen? 
In  order  to  give  a  (pre I i m i nary)  answer  to  these  questions  we  shall  lake 
a  closer  look  at  spherical  multipole  layers  in  the  next  subsection. 

7.3.  Multipole  layers  on  the  sphere. 


Assume  a  sphere  of  radius  R  =  I  and  let  r  >  R  =  I.  Starting  from  the 
well  known  representation  for  the  function 


|  with  l  =  -\J~R  +  r  -  2  R  r  cos  | 
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Bfl 


namely. 


I  I 


l 


I  Pn(cos  f)  {^}n  ;  r  >  R 


(7.6) 


n  =  8 


differentiating  N  times  with  respect  to  R  and  putting  afterwards  R  =  1, 
one  gets 


r__  a 

dR 


f  =  I  P  (cos  |)  n(n-1 ) . . . (n-N+l ) 

1  n=N  n  r  rn+1 


(7.7) 


Here  Pn(t)  are  the  familiar  Legendre  polynomials,  normalized  in  a  way 
that  Pn(D  =  I  -  The  series  (7.7)  converges  for  r  >  1.  Letting  i — “I,  we 
obtain  in  the  Ijmil  a  sequence  which  converges  distr ibut ional ly .  The 
potential  generated  by  the  multipole  layer  with  strength  function  fi  is 


V<rf>  =  I 


\ 


r  3rn 


ix(tj)  drc-r,) 


(7.8) 


>  R  =  1 


Here  £,  i)  are  unit  vectors,  T  is  the  surface  of  the  unit  sphere.  We 
temporarily  fix  r  and  view  V(r£)  as  a  function  of  f  only.  Then  (7.8) 
represents  an  isotropic  operator.  According  to  the  theory  outlined  in 
Meissl  (1971),  chapter  3,  we  know  that  its  eigenfunctions  are  the 
spherical  harmonics  Hnm.  We  can  calculate  the  eigenvalues  by  the 
Funk-Hecke  formula  obtaining 


=  0 


;  n  <  N 


-2^r  ;  "  *  N  “  0 

4* 


(7.9) 


2n  +  I 


7  n(n~! ) . . . (n-N+l ) 


n+  I 


;  n  5=  N  >  8 


at 
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Hence  if  (xCJ)  is  represented  as 

ftp  ■  e„.  Hn.c?}  " 

n  =  0  m  -  -  n 

where  we  now  take  the  sur  face~spher  i  ca  I  harmonics  H^Ci))  as  fully 
normal ized,  i .e. 


a tfr  W  m  1 


we  get  the  representation  for  V(r£)  as 


»  +  n  H  (  £  ) 

VCrf>  -  l  »n  ?  > 

n  -  N  m  —  r* 


im  r  "  *  1 


(7.12) 


We  see  that  VC r?)  =  O^},  as  announced  earlier.  Further  we  note  that 
the  operator  is  singular.  The  coefficients  finm  for  n  <  N  do  not 
contribute  to  V(r£). 

If  we  like  to  continue  V(r£)  downward  to  r  =  R  -  1/  we  must  impose 
restrictions  upon  the  jinm .  The  following  series  must  converge 


I 

n  =  N 


+  n 

I 


m  =  -  n 


<  « 


(7.13) 


to  ensure  that  V(f)  is  a  member  of  HP,  the  Hilbert  space  of  squared 
integrabl e  functions  defined  on  T.  The  operator  transforming  /x(£)  into 
VCJf )  has  eigenvalues 

%  =  n(n-1 )  ...  (n-N+l ) 

n  2n+ 1 


It  amplifies  the  high  frequences  more  than  the  lower  ones.  It  may  be 
viewed  as  an  operator  from  the  Hilbert  space  H^  N  into  H^;  k>0. 


comprises  those  functions 


The  space  Hp^ 

f(J)  =  l  l  c  H(J)  C7.I5) 

j  n  nm  nm  ^ 

n  =  o  m  =  -  n 

for  which  the  following  series  converges. 

I  T  (n'c  )’  <  »  a  l6) 

_ «  nm 

n  =  U  hi  =  -  n 

In  a  least  squares  approach  we  like  to  have  at  r  s  I  a  potential 
V(£)  e  Hp25 .  Hence  |i(£)  e  HpN+l)  appears  appropriate.  Therefore  the 
trial  functions  should  be  CN  across  element  boundaries.  In  a  Ritz-type 
approach  V(£)  e  Hp1^,  fi(£)  6  HpN5  and  CN-1  continuity  across  element 
boundaries  would  be  sufficient. 

Remark^  One  could  avoid  the  smoothness  requirements  on  fi(£)  by  the 
Bjerhammai — sphere  approach.  A  sphere  of  radius  |-e  is  placed 
concentrically  into  the  unit  sphere.  The  density  \l  is  assumed  on  this 
sphere.  We  thus  consider  an  operator 

|lC(l-c)J)  —  vcp  (7.17) 

whose  eigenvalues  are 

=  2n+T  nCn'° 

However,  in  the  present  context  we  do  not  propose  to  shove  difficulties 
under  the  carpet  by  burying  the  mullipole  layer  underneath  the  surface 
of  computation 

Suppose  now,  for  the  moment,  that  in  (7.8)  VCrjp  is  specified  and 
that  fi(‘T})  is  sought.  If  the  equations  are  discretized  in  some  way,  we 


n  -  N 


...  Cn-N+1  XI  -e) 

n  ?  N 


(7.18) 
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expect  them  to  be  very  ill  conditioned  due  to  the  differences  in  size  of 
the  eigenvalues  X^  .  While  the  high  frequent  components  of  fiCi))  appear  to 
be  well  determined,  the  low  frequent  components  are  not. 

In  physical  geodesy  we  do  not  start  from  a  known  VCJ).  Instead  we 
have  a  discrete  set  of  measurements  corresponding  to  functionals  of 
V(£).  The  situation  is  more  difficult  but  nevertheless  similar.  The 
linear  system  leading  to  will  be  ill  conditioned. 

One  may  adopt  the  viewpoint  that  ill  conditioning  is  equivalent  to 
an  inproper  problem  formulation,  and  one  may  simply  abandon  the  outlined 
approach.  On  the  other  hand,  if  insight  teaches  one  that  the  effects  of 
the  ill-conditioning  on  the  final  result  will  be  negligible,  one  may  find 
a  way  around  the  numerical  difficulties.  Of  course,  there  is  no  hope 
that  insight  can  be  replaced  by  calculus. 

At  the  present  time,  no  definite  answer  can  be  given.  More  lime 
is  needed  to  think  about  the  method  and  to  conduct  numerical 
exper imenls. 
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Appendix  A 

Computat i ona 1  effort  associated  with  the  partial  reduction  of  o  profiled 
system  of  normal  equal  ions. 

The  material  of  this  section  is  rather  unsophisticated  in  nature. 

However  even  trivial  considerations  may  be  cumbersome  if  their 
verification  is  left  to  the  reader.  Therefore  we  quickly  state  some 
formulas  on  the  number  of  operations  needed  to  solve  a  structured  large 
system  of  positive  definite  equations.  A  more  systematic  introduction  to 
this  problem  area  is  found  in  Meissl  Cl 980  ),  chapter  6. 

We  introduce  the  well  known  concept  of  the  profile  of  a  symmetric 
matrix  A  =  (a^).  It  comprises  all  elements  aLj  such  that  Cl)  i  s  j,  and 
C2)  there  exists  an  element  akj  +  0,  k$i.  Referring  to  figure  A.  I  we 
introduce  the  profile  function  pCx)  indicated  by  a  heavy  solid  line.  The 
heavy  line  should  actually  be  a  step  function  because  we  have  a  discrete 
number  of  equations.  However  we  smear  out  the  discontinuities.  This  is 
legitimate  if  we  deal  with  a  large  system. 


Figure  A. I.  Partial  reduction  of  a  profiled  symmetric 
matrix.  Definition  of  pCx),  dCx). 
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The  portion  below  the  main  diagonal  is  not  shown  in  figure  A.!.  Due  to 
symmetry,  there  is  no  need  to  store  these  coefficients  in  the  computer. 

The  upper  portion  of  the  system  comprises  i  equations.  They  are 
called  “interior".  The  equations  below  labeled  j  are  called  "junction 
equations’.  During  partial  reduction,  the  interior  equations  are 
eliminated.  One  obtains  a  system  of  j  partially  reduced  junction 
equat i ons . 

If  the  matrix  is  split  as 


A  * 


CA.1) 


then  the  matrix  of  the  partially  reduced  set  is  given  by 


A 


(p) 

ii 


(A. 2) 


A  direct  elimination  procedure  is  used  such  as  Gauss,  Cholesky  or  one  of 
the  variants.  The  partial  triangular  decomposition  phase  is  by  far  most 
time  consuming.  During  this  phase  zero  coefficients  are  enforced  below 
the  main  diagonal  positions  of  the  interior  equations.  Focus  attention 
on  a  row  of  coefficients  whose  diagonal  position  is  implied  by  column  x. 
Call  this  the  row  x.  Let  y  be  a  column  to  the  right  of  x.  From  the 
coefficient  in  row  x  and  column  y  as  many  multiples  of  coefficients 
located  above  it  are  subtracted  as  the  following  expression  indicates: 

riKY)  *  Min  {  Mar(  p(x)~  d(x)f  o)t  H av  (  piy) -  dCx), o)}  (A. 3) 


Thereby  dCx)  is  the  function  implied  by  the  heavy  broken  line  in  figure 
A.I.  Hence  partial  triangular  decomposition  requires 

F  ■■  J cU  j P  (*,  yj  oijr 

o  n 


CA.4) 
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steps.  One  step  comprises  one  multiplication  and  one  addition.  There 
are  also  divisions  and  square  roots  Cin  case  of  Cholesky),  but  their 
number  is  negligible.  Also  the  number  of  operational  steps  involving 
the  right  hand  side,  as  well  as  the  number  of  steps  arising  during  the 
later  back-substitution  phase  are  negligible. 

Take  a  fully  occupied  system.  We  find 


L  L+J  Lt J  .  l 

r  *  J  d.x JycLy  +  J  -fUy  s  ~r  + 


LJ  (L+j) 


(A. 5) 


This  formula  is  frequently  used  in  chapter  4.  For  j  =  0,  i.e.  full 
reduction,  we  obtain 


r 


(A. 6) 
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